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C D , We review the topic of quantized vortices in multicomponent Bose-Einstein condensates 

of dilute atomic gases, with an emphasis on that in two-component condensates. First, 
' we review the fundamental structure, stability and dynamics of a single vortex state in a 

slowly rotating two-component condensates. To understand recent experimental results, 
we use the coupled Gross-Pitaevskii equations and the generalized nonlinear sigma model, 
i-^j ■ An axisymmetric vortex state, which was observed by the JILA group, can be regarded 

(-H ' as a topologically trivial skyrmion in the pseudospin representation. The internal, coher- 

Q | ent coupling between the two components breaks the axisymmetry of the vortex state, 

, resulting in a stable vortex molecule (a meron pair). We also mention unconventional 

vortex states and monopole excitations in a spin-1 Bose-Einstein condensate. Next, we 
discuss a rich variety of vortex states realized in rapidly rotating two-component Bose- 
k>( ' Einstein condensates. We introduce a phase diagram with axes of rotation frequency 

, and the intercomponent coupling strength. This phase diagram reveals unconventional 

^ i vortex states such as a square lattice, a double-core lattice, vortex stripes and vortex 

sheets, all of which are in an experimentally accessible parameter regime. The coherent 
coupling leads to an effective attractive interaction between two components, providing 
not only a promising candidate to tune the intercomponent interaction to study the 
rich vortex phases but also a new regime to explore vortex states consisting of vortex 
molecules characterized by anisotropic vorticity. A recent experiment by the JILA group 
vindicated the formation of a square vortex lattice in this system. 

Keywords: Bose-Einstein condensation; superfluidity; quantized vortex; vortex lattice; 
binary system; rotation; spinor condensate; spin texture; relative phase 
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1. Introduction 

Due to their being a common feature of fluid motion, vortices appear in diverse fields 
of science. In a superfluid, that is, a fluid without viscosity, the flow has quantized 
circulation which originates in quantized vortices 1 . A quantized vortex is commonly 
viewed as a topological defect of the order parameter in Bose-Einstein condensation. 
Since the gradient of the phase of the order parameter is proportional to the super- 
fluid velocity, quantized vortices play a crucial role in understanding superfluidity. 
Moreover, the structure of the order parameters in multicomponcnt Bose-Einstein 
condensates, hereafter BECs, are analogous to other studies including superfluid 
3 He 2 , unconventional superconductors 3 , quantum Hall systems 4 , nonlinear optics 5 , 
nuclear physics 6 , and cosmology 7 . Therefore, the investigation of unconventional 
topological defects in such multicomponent BECs is of great importance and of 
broad interest. 

While quantized vortices have been extensively studied in the field of superfluid 
helium 1 , they have seen a remarkable resurgence of interest since the realization 
of Bose-Einstein condensation in dilute alkali-atomic gases 8 ' 9 . This Bose-condensed 
system has some distinct advantages for the study of quantized vortices. First, the 
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BECs of alkali-atomic gases are so dilute that the interatomic interaction can be 
accurately parametrized in terms of the s-wavc scattering length. Experimentally, 
the diluteness leads to a relatively large healing length that characterizes the vor- 
tex core size, thus enabling us to directly visualize the quantized vortices by optical 
techniques. Also, because the gas is dilute, the Gross-Pitaevskii equation gives an ac- 
curate description of static and dynamic properties of the atomic condensate within 
a mean-field approximation. This situation differs from superfluid 4 He in which the 
high density and the strong interaction complicate the theoretical treatments. Sec- 
ond, rapidly developing techniques of atomic, molecular and optical physics allow 
us to carry out quantum-state engineering of the condensates in a well-controlled 
manner, thereby opening a new direction for the study of coherent quantum phe- 
nomena. For example, since alkali atoms have internal degrees of freedom attributed 
to the hyperfine spin, an external field can be used to induce conversions between 
spin states and thus it can be used to control the spatial variation of the wave 
functions. Moreover, a field-induced Fcshbach resonance allows one to tune of the 
atom-atom interaction and create of a quasi-bound molecular state 10 . 

The first experimental detection of a vortex in a dilute alkali-atomic gas BEC 
was made by Matthews et al. in 1999 11 . Their study involves condensates of 87 Rb 
atoms residing in two hyperfine states. The generation of a vortex was achieved by 
a phase imprinting technique proposed by Williams and Holland 12 , in which the 
interconversion between two components was controlled spatially and temporally 
by an external coupling field. After turning off the coupling, the vortex state is 
formed by one circulating component surrounding a nonrotating core of the other 
component. Leanhardt et al. 13 used "topological phase imprinting" 14 to create a 
vortex in a spin-1 BEC. These two methods use internal degrees of freedom of the 
condensate to create vortices, but a simple way to create vortices is to stir a con- 
densate mechanically with a rotating "bucket". Following the latter method, the 
ENS group observed the formation of a single vortex and multiple vortices in a 
single-component elongated condensate 15 . The condensate was trapped in a static 
axisymmetric magnetic trap which was deformed by a nonaxisymmetric attractive 
dipole potential created with stirring laser beams. This combined potential pro- 
duces a cigar-shaped harmonic trap with a slightly anisotropic transverse profile. 
By rotating the orientation of the transverse anisotropy at frequency fi, they ob- 
served the dynamic formation of a vortex above a certain critical value of 0. A 
lattice involving more vortices appeared when O was further increased. In contrast 
to indirect visualization methods used in superfluid helium systems, the quantized 
vortices were directly visualized as hollows in the transverse density profile of the 
time-of-flight image. Later, groups at MIT 16 , JILA 17 , and Oxford 18 observed vortex 
lattices with the first two groups succeeding in creating 100 or more vortices. 

Theoretical studies mainly used the mean-field Gross-Pitaevskii theory to de- 
scribe the main features of the vortex states (see Ref. 19 for a review), and several 
predictions have been shown to agree with experiments. Some of the important 
problems are concerned with the equilibrium properties of a single vortex, includ- 
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ing its structures and dynamics 20 - 24 , the critical frequency and the surface- mode 
instability for vortex nucleation 25-33 , and the nonlinear dynamics of vortex lat- 
tice formation 34-39 . Finally, experiments on a rapidly rotating BEC 40-44 have mo- 
tived theoretical studies of collective oscillations of a vortex lattice 45-50 , new vortex 
phases in an anharmonic potential 51-59 , and a strongly-correlated fractal quantum 
Hall phase via quantum melting of a vortex lattice 60-67 . 

Another important issue in vortex physics is to elucidate the vortex phase in 
multicomponcnt BECs. Multicomponent order parameters allow the formation of 
various unconventional topological defects with complex properties that arise from 
interactions between different order-parameter components 68 . Since it is possible 
to load and cool atoms in more than one hyperfine spin state or more than one 
atomic element in the same trap, multicomponent condensates can be realized 
experimentally 69-77 . Such systems offer an ideal testing ground for the study of 
such topological defects; moreover the microscopic characterization of these sys- 
tems should in principle be possible because this system is free from impurities and 
well controlled by optical techniques. 

In this paper we review recent studies on quantized vortices in multicomponent 
BECs, especially focusing on those in two-component BECs. To the best of our 
knowledge, there are no reviews of two-component BECs despite the large number 
of studies of this system 78-96 . Also, while a number of interesting phenomena have 
been reported on vortices of a spin-1 BEC, we mention this system only briefly and 
would like to refer the reader to references 97-117 for more details. 

This paper is organized as follows. After a brief introduction of quantized vor- 
tices in a single-component BEC in Sec. 2, we formulate in Sec. 3 a mean-field 
theory that describes two-component atomic-gas BECs and give a brief review on 
the ground state structure without vortices. In Sec. 4, after describing experiments 
on the creation and observation of vortices in two-component BECs 11 , we discuss the 
static structure and dynamic properties of the single- vortex state in two-component 
BECs in a slowly rotating regime. By introducing the "pseudospin" that describes 
the two-component BECs, we discuss the properties of topologically "trivial" and 
"nontrivial" spin textures. In particular, we predict a vortex molecule in a two- 
component BEC with an internal coherent coupling. We also mention briefly un- 
conventional vortex states in a spin-1 BEC. In Sec. 5, we consider the vortex states 
of two-component BECs under a rapid rotation. Combination of an analytical ap- 
proach and numerical simulations reveals a rich variety of vortex states, which is 
not found in a single-component BEC. We also discuss the effect of the internal 
coherent coupling on the structure of vortex states. Finally, we conclude this review 
article in Sec. 6. 

2. What is a quantized vortex? 

Below a critical temperature in an ideal Bose gas, a finite fraction of the particles 
occupies the same single- particle ground state and forms a BEC. When the particles 
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interact, the single-particle density matrix n^(r,r') = (W(r)^>(r)) is still defined, 
where W is the boson creation (anihilation) operator. When a single-particle 
state is occupied by a macroscopic number of particles, n^(r, r') has a macroscopic 
eigenvalue. We may define a condensate wave function ^(r, t) as the corresponding 
eigenstates 118 . In the absence of interactions this wave function reduces to that 
of the macroscopically occupied single particle state. When a BEC is held in an 
external potential V(r), the dynamics of ^(r, t) is described by the Gross-Pitaevskii 
(GP) equation 8 ' 9 

= (-|!-V 2 + V(r) + g\Hf(r, t)| 2 ) *(r, t), (1) 

where g = A-Kh 2 a/m represents the strength of interaction characterized by the 
s-wave scattering length a and particle mass to. The condensate wave function is 
normalized by the total particle number N as j d 3 r\^>\ 2 — N. Writing \t(r, i) = 
\J p(r, t) exp(i#(r, t)), the absolute squared amplitude p(r,t) = |^(r, t)\ 2 gives the 
condensate density and the gradient of the phase 9{v, t) gives the superfluid velocity 
v s (r,i) = (h/m)V6(r,t). As a result, the circulation of v s around a closed path C 
in the fluid is quantized as 

/ ds ■ v s = — / ds ■ V0 = uk (n = 0,±1,±2, . . .) (2) 
Jc m Jc 

with the quantum of circulation k — h/m. Such a vortex with a quantized circulation 
is called a "quantized vortex" . 

Let us review the fundamental properties of quantized vortices. Consider a single 
straight vortex line along the symmetry axis of a cylindrical vessel that can rotate 
at frequency O = ilz. In this case, the velocity field of the vortex line is 

v S = ^ (3) 

where r is the radial distance from the line and 9 is the azimuthal angle. The 
corresponding "vorticity" is 

V x v s = K<5 (2) (r ± )z. (4) 

There is a singularity on the vortex line (r = 0), where v s diverges. Near r = the 
superfluid density is significantly supressed, forming the vortex core. In the case of 
superfluid helum, the size of the vortex core r c is of the order of the interatomic 
distance, while in a dilute gas this size equals the healing length r c ~ £, where 

(5) 



V 2m 9P 

The angular momentum L z of the fluid due to the vortex line is 

f , -kLR 2 
L z = pv s rdr = np^—, (6) 

where L is the height of the vessel, R is its radius, and p is the value of the superfluid 
density far from the core. In deriving the last equality we assumed r c <C R- 
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The energy associated with the vortex line is dominated by the kinetic energy 
term. A straight calculation yields 

E v = J \ P v 2 s dv^lL j\^rdr = n 2 L^\n{^ . (7) 

Since E v oc n 2 , vortices with n > 1 are energetically unfavorable. That is, the energy 
cost to create one n = 2 vortex is higher than that to create two n = 1 vortices, even 
though they give the same angular momentum for R 3> 1 (L z oc n). Therefore, a 
stable quantized vortex always has n = 1; the disintegration of a multiply quantized 
vortex has been intensely studied 119-122 and it was observed experimentally by Shin 
et al. in an atomic-gas BEC 123 . 

Using Eqs. (6) and (7) one can easily find the critical frequency of rotation f2 c 
for the existence of an energetically-stable vortex line. In a frame rotating with 17, 
the free energy is given by F = E — ilL z and the presence of a vortex increases E 
by E v . For a vortex with n = 1, the critical frequency is obtained by imposing that 
the change in the energy QL Z due to the creation of the vortex be equal to E v as 

For 17 > 17 c the free energy of the system with one vortex is lower than that 
without a vortex. When the rotation frequency is significantly higher than 17 c , more 
vortices will appear in the form of a triangular lattice of vortices. For very large 
17 the rotation of the superfluid will then look similar to the one of a rigid body 
characterized by V x v s = 217. Using Eq. (4), one finds that the average vorticity 
per unit area is given by V x v s = nn v z, where h v is the number of vortices per 
unit area, so that the density of vortices is related to the rotation frequency 17 as 

217 

n v = — . (9) 

K 

This relation, often referred to as Feynman's rule 124 , can be used to determine the 
maximum possible number of vortices in a given area as a function of 17. 

The above formalism is extended to quantized vortices in a trapped BEC by 
including the explicit form of the trapping potential V(r). This potential leads to 
an inhomogeneous condensate density, which plays an important role in the stability 
and dynamical behavior of the vortices. As an example, consider the structure of 
a single vortex in a condensate trapped by the axisymmetric harmonic potential 
V(r,z) — muj' 2 _(r' 2 + X 2 z 2 )/2 with the aspect ratio A = uj z /oj±. The equilibrium 
state is determined by the time-independent GP equation 

-|^V 2 + V(v) + .g|* (r)| 2 ) *„(r) = M*o(r), (10) 

where the time dependence of *o is fixed by the chemical potential yu as 
*(r, t) = * (r)e- i ' i */'\ The value of /i is determined by the normalization condition 
/ | 'Jo ( r ) 1 2 = N. The condensate wave function with a singly-quantized vortex line 
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located along the z-axis takes the form \&o( r ) = f(r, z)e t9 with the cylindrical coor- 
dinate (r, 0, z), where / is a real function that is related to the condensate density 
as po = f 2 . Then, Eq. (10) becomes 

/r ( <>■' i d <■>■■■ \ n 2 . 

/ = m/- (11) 



2m \dr 2 r dr dz 2 ) 2mr 2 ^ ^ 

Here, the centrifugal term H 2 /2mr 2 arises from the azimuthal motion of the con- 
densate which gives rise to a kinetic energy density pomv 2 /2 = p h 2 /2mr 2 . 

The trapping potential defines the characteristic length scale &ho = \Jh/muj^. 
Without an interparticle interaction (g = 0), the scale &ho determines the spatial 
extent of the wave function. An increase in the repulsive atomic interaction expands 
the condensate, so that the kinetic energy associated with the density variation can 
be made negligible compared with the trap energy and the interaction energy The 
ratio of the interaction energy to the harmonic-oscillator energy gpo/fiu!± ~ Na/bh 
serves to quantify the effect of the interaction. In the limit Na/bho S> 1, which 
is applicable to current experiments on trapped BECs, one can omit the terms 
involving derivatives with respect to r and z in Eq. (11). This approximation is 
called the Thomas- Fermi approximation 125 , and the density profile for the vortex 
state is then given by 



. . If mbj 2 , . 9 , 99n h 2 \ ^ ( mid 2 , , 9 ~ 9 , h 2 
Po(r, z) = -[p--^ { r 2 + X 2 . 2 ) -—} @ ^-—L { r 2 + X z ) - — 



roWi 7-2 z " €(0)2 Wi r2 z " m '\ (12) 
-p°W\i- 1 ^-zi-—)o{i-xr- — J, (12) 



where Q(x) is a step function and po(0) is the density at the center of the vortex- 
free Thomas-Fermi profile. We define the Thomas-Fermi radius R± = yj2p/muj 2 _ 
and R z = y / 2p/mu; 2 and the healing length £(0) = y / h 2 /2mgp a (0). Due to the 
centrifugal term £(0) 2 /r 2 , the condensate density vanishes at the center out to a 
distance of order £, whereas the density in the outer region has the form of an 
upward-oriented parabola. The Thomas-Fermi approximation is valid in the region 
£ < r < R±. The asymptotic form of p (r) for r < £ is 126 p (r) = p o (0){r /£{0)) 2 
whose interpolation to Eq. (12) gives 

r 2 \ ( r 2 z 2 \ „ / r 2 z 2 



p^z) * po(o) [i - w± - W J e\i - ^ - M ) . (i3) 

This form agrees fairly well with the profile of the single vortex state obtained by 
the numerical calculation of the GP equation 126 . 

For a detailed account of the single vortex state, we refer to a review article by 
Fetter and Svidzinsky 19 (see also Chap. 9 in Ref.8 and Chap. 14 in Rcf.9). 
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3. Formulation of two-component Bose-Einstein condensates 
3.1. Coupled Gross-Pitaevskii equations 

We now consider two-component BECs as described by the condensate wave func- 
tions for the two hyperfine spin states and xf^ 11 ' 69 ' 70 - These wave functions can 
also describe a system of two atomic species 77 ' 127,128 or two isotopes of an atomic 
element. The theory for the latter case is described in Refs. 129 and 130, and ex- 
periments are described in Ref. 131. In addition, the condensates are confined by 
trapping potentials Vi(r) (i = 1,2), which is assumed to rotate at a rotation fre- 
quency f2 about the z axis as fl = fiz. Viewed from the rotating frame, the GP 
energy functional of the two-component BEC is 8 ' 9 



i=l,2 ^ m * 



(14) 



Here, rm is the mass of the i-th atom and VtL z = ihil(yd x — xd y ) is a centrifugal 
term. The coefficients (i = 1,2) and gi2 represent the atom-atom interaction 
between atoms iin the same hyperfine state, and atoms in different hyperfine states, 
respectively. They are expressed in terms of the corresponding s-wave scattering 
lengths ai, ci2, and a\i as 

47rfi 2 a l 

ft = , (15 

m,i 

2TrH 2 a 12 . , 

5i2 = , (16) 

mi2 

where — "m-i 1 + is the reduced mass. 

When the trapping potential Vi{r) is produced by a magnetic field, each com- 
ponent has its own potential energy because the potential depends on a magnetic 
moment of the atom. For example, the trapping potential for atoms with the hy- 
perfine states \Fi,Mpi) is expressed as a function of the magnetic field \B(r)\ as 127 



Vi{r) = u B GiM Fi \B{r)\ ~ n B GiM Ft 



B + -(K x x 2 + K y y 2 +K z z 2 ) 



= V° + ^mi(wf x x 2 + u; 2 y y 2 + lj 2 z z 2 ), (17) 

where u B is the Bohr magneton, Gi is the g-factor of the i-th atom, V® = 
/j B GiM Fi B , and uJik (k = x,y,z) is the trapping frequency satisfying the rela- 
tion 

m t ujj k =u B GiM Fi K k . (18) 

The potential minima V t ° for the BECs with different values of GiM Fi do not 
coincide. This positional shift is also influenced by a gravitational effect, by the dif- 
ference of nuclear magnetic moments, and by nonlincarity in the Zecman shift 69 ' 70 . 
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If the potential is created by a pure optical laser 71 ' 73 ' 74 ' 75 , all atoms share the same 
harmonic potential. 

The time-dependent coupled GP equations for two-component BECs can be 
obtained by using a variational procedure ihd^i = SE/S^* as 

city h 2 

ih—± = -—V 2 *! + Ki(r)*i + ff i|*i| 2 *i + 3i 2 |* 2 | 2 *i - QL X 9 U (19) 
at 2mi 

Oty h 2 

= -— V 2 * 2 + y 2 (r)* 2 + <fc|* 2 | 2 *2 + 0i 2 |*i| 2 *2 - nL»*2. (20) 
ot 2m 2 

Substituting *j(r, i) = \I/j(r)e _l ' ii */'' yields the time-independent coupled GP equa- 
tions 

h 2 

V 2 *i + ^i(r)*i + 5 i|*i| 2 *i + ff i 2 |1'2| 2 *i - = A*i*i, (21) 

V 2 * 2 + V 2 (r)^ 2 + .92|*2| 2 *2 + fl r 12|*l| 2 *2 ~ fi£**2 = ^2*2, (22) 



2mi 
ft 2 



2m 2 

where we have introduced the Lagrange multiplier fa which represents the chemical 
potential and is determined so as to satisfy the conservation of particle number 
Ni = J d 3 r\^i\ 2 . In the two-component BECs, each component interacts with the 
other through the intercomponent mean- field coupling cx 312, which yields structures 
and dynamics not found in a single-component BEC. 



3.2. Ground state structure 



It is instructive to start with the ground state of two-component BECs without ro- 
tation (O = 0). The equilibrium solutions of the coupled GP equations (21) and (22) 
show a rich variety of stable structures 78 ' 130 ' 127,132 " 142 , depending on the various 
parameters of the system. In particular, the intercomponent interaction g i2 plays an 
important role in determining the structure of the ground state; when the intercom- 
ponent interaction is strongly repulsive, the two components phase separate 136 ' 137 . 
The qualitative features of the ground state for a trapped two-component BEC can 
be understood by using the Thomas-Fermi approximation 78 ' 127 ' 140 in which one ne- 
glects the quantum pressure terms (h 2 /2rrii^/pi)V 2 ^fpl with the condensate density 
pi = \^i\ 2 in Eqs. (21) and (22). First consider the simple case G\M Fl = G 2 Mp 2 
in Eq. (18) (i.e., m\U)\ k — m 2 uj' 2 k )- The equilibrium density distributions p- ^ in the 



spatial region where the two components coexist (p\ 



(0) 



^ and p ( 2 ] 



(o) 
Pi 



92P1 - 912P2 - 92V1 + g 12 V 2 



1 



(0) 91P2 
P2 = 



3i52 
912H1 



9u 



9iV 2 +gi 2 Vi 



3i r i2 
1 



Mi 



.912 
92 



-M2 



^ 0) are 

r 2 Vi(r)L 



3152 - 9{ 2 92^X2 

where, we defined the following dimcnsionless parameters 



M2 - —Mi - riV2(r) 
9i 



12 



1 



mi + m 2 a{ 2 
4m 12 O1O2' 



(23) 
(24) 

(25) 
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TOl ai2 
1 1 = 1 _ T, ' 

2mi2 ai 

_ m 2 ai2 

1 2 = 1 - 7: • 

2mi2 a 2 

For the two components to coexist, the right-hand sides of Eqs. (23) and 
be positive. Then, the parameter Ti 2 must be positive, that is, 

3i52 - 3i2 > 0. 

The coexisting region is given by the region delimited by the surface defined as the 
solution to the equation yn — gi2Pj/gj — TjVi(r) = 127 . Otherwise, that is when 

3i32 - 3r2 < 0, (29) 

there is no coexisting region, which implies that the two components phase-separate. 

For Ti 2 > 0, there are two possible density profiles that depend on the signs of 
the parameters Ti and r 2 . These signs determine the curvature of the density profile 
at the center of the trapping potential because Ti and T2 are the coefficients of the 
x 2 , y 2 and z 2 terms in Eqs. (23) and (24). For Ti > and T 2 > 0, both density 
profiles are written by upward-convex parabolas in the coexisting region. For Ti > 
and F 2 < (ri < and T 2 > 0), one density profile is an upward-oriented parabola, 
but the other is a usual parabola. In the peripheral region in which the density of one 
component vanishes, the other component follows the single-component Thomas- 
Fermi profile p\°' = (fii — Vi(r))/gi. This profile is connected to the density profile 
of Eqs. (23) and (24) by adjusting the chemical potential pi (i = 1 or 2) subject to 
the normalization condition J d 3 rpf*^ — Ni 78 . 

For Ti 2 < the two components phase-separate as the coexisting regions do 
not exist. In the region in which only one density is nonvanishing (p\ ^ ^ and 
p^ = 0, for i 7^ k), the density profile is given by the single-component Thomas- 
Fermi profile = (pi — Vi(r))/gi. Since the Thomas-Fermi approximation fails 
to describe the domain boundary region, where the quantum pressure term cannot 
be neglected, numerical analysis of Eqs. (21) and (22) is necessary to determine the 
explicit density profile in this regime 135 . The numerical calculation shows that there 
are two characteristic configurations of the phase-separated two-component BEC in 
the trapping potential. One configuration consists of a core of one component at the 
center of the trap and a shell of the other around the one component, a structure 
that preserves the spatial symmetry of the trapping potential 78 ' 135 . The other con- 
figuration breaks the spatial symmetry by displacing the center of each condensate 
from the center of the trapping potential 134 ' 138 . Whether the ground state takes the 
former or the latter configuration depends not only on the intercomponent interac- 
tion but also on the particle number, intracomponent interaction, and the shape of 
the trapping potential 139,141 ' 142 . 

Compared to single-component BECs, two-component BECs exhibit a rich va- 
riety of ground state structures. These structures and their stability depend upon 
various parameters, particularly on the condition of the phase separation given by 



(26) 
(27) 
(24) must 

(28) 
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Eq. (29). Two-component mixtures of BECs with different spin states have been 
created in the laboratory for the systems of 87 Rb (Refs. 69, 70) and 23 Na (Ref. 72). 
For the case of 87 Rb, the nearly-equal scattering lengths (ai = 5.67 nm, a 2 = 5.34 
nm, and a\i — 5.50 nm) and the small relative displacement of the trapping min- 
ima lead to a weak separation between the two components. In contrast, for the 
23 Na-mixture, the value a\ — 2.75 nm, a 2 = 2.65 nm, and a\i = 2.75 nm causes a 
spin domain to form 71 ' 72 . Current experiments are aiming to create two-component 
BECs with two atomic species 128 ' 143 ; a mixture of the 41 K and 87 Rb BEC has been 
observed 77 . In this case, the intercomponent scattering length a 12 may take val- 
ues different from the intracomponent ones a± or a 2 . Furthermore, we should note 
that an additional possibility to change the value of a\ 2 by the use of the Feshbach 

129 144— 147 

resonance J ' '. 

4. Vortex states in slowly rotating two-component BECs 
4.1. Creation of a quantized vortex 

Although most theoretical studies on quantized vortices have focused on single- 
condensate systems, the first observation of a quantized vortex in a gaseous BEC 
was achieved in a two-component system 11 ' 12 . In the JILA experiment 11 , the vortex 
was created in a two-component BEC consisting of 87 Rb atoms with the hypcrfinc 
spin states \F = 1, M F = -1) = |1) and \F = 2, M F = 1) = |2). Since the scattering 
lengths between atoms of |1) and |1), |2) and |2), and |1) and |2) are all different, the 
two states are not equivalent, and the two-component condensate is characterized 
by the two-component order parameter. 

The advantage of using the |1) and |2) states of 87 Rb atoms to perform quan- 
tum state engineering is that, first, their magnetic moments are nearly equal, so 
that they can be confined simultaneously in almost identical magnetic potentials. 
This is advantageous to investigate the intrinsic interaction phenomena between two 
condensates. Second, nearly equal s-wave scattering lengths for two 87 Rb atoms in 
|1) and 1 2) result in the fortuitously small inelastic spin-exchange rate 148 , which 
provides a stable binary system with a long lifetime. In contrast, the inelastic rate 
is very large in 23 Na atoms. Third, we can use a microwave rf-pulse to change the 
populations of atoms from the |1) state to the |2) state 11 ' 149,150 . If a microwave is 
applied continuously, this system is referred to as an "internal" Josephson effect 
because the coexisting components are coupled through their internal degrees of 
freedom 151-154 . Thus, this system offers new possibilities not only for quantum- 
mechanical state engineering, but also for exploring new vortex structures in coher- 
ently coupled multicomponent BECs as described in Sec. 4. 3. 

In experiments at JILA experiments , atoms are initially trapped in one state, 
say, the |1) state, and then cooled below the Bose-Einstein transition temperature. 
When the atoms in the |1) state forms a condensate, a two-photon microwave field 
is applied, inducing transitions between the |1) state and the |2) state. As a result, 
the atoms undergo coherent oscillations between the two hyperfme levels with an 
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effective Rabi frequency w^ ff = \J + A 2 , where the Rabi frequency wr gives 
the rate at which population would oscillate between the two states when A = 
and depends on the microwave power and the atom- field coupling constant. The 
detuning parameter A denotes the difference between the frequency of the coupling 
rf field and the frequency difference between the two internal atomic states. 

The effect of internal coupling can be formulated by introducing the energy 
density 

E c = -hLu R (V* 2 y ie - lAt + *J* 2 e lAt ), (30) 

which allows atoms to change their internal state coherently. Since the microwave 
coupling fields are time dependent, it is useful to use the frame in which the field is 
constant (a frame rotating with the laser field) . Then, the coupled time-dependent 
GP equations read 

ih ^ir = f-^r^ + V i+5il ,I, i| 2 + ffi2l*2| 2 + ^)*i-^'R*2, (31) 
at \ Ira I ) 

lh ^f = (-^+ y 2+32|*2| 2 +5l2|*l| 2 -^)*2-^R*l. (32) 

The last terms in Eqs. (31) and (32) cause interconversions between the components. 
For a homogenous system without Vj but having uniform phases of the condensate 
wave functions for both components, the interconversion takes place at the same 
rate everywhere. However, the time variation of the spatially inhomogeous potential 
Vi changes the way the interconversion is made, which is a key ingredient of the 
phase imprinting method for vortex creation 12 . 

To proceed, let us assume that the external potential consists of a radi- 
ally isotropic harmonic potential and a time-dependent perturbation V\(r, t) — 
muJ]j 2 /2 + H x and V 2 (r, t) = muj 2 L r 2 /2 - H 1 with 

Hi = K [/(r) cos(ft'i) + g(r) sin(fi'i)], (33) 

where k characterizes the strength of the perturbation and /(r) and g(r) describe 
the spatial dependence of the perturbation. The form of H\ can be determined from 
the symmetry of the quantum state to be prepared. To create a vortex state with 
one unit of angular momentum, one can take n — muj\r ai /(r) = x, and g(r) = y 
in Cartesian coordinates 12 . This form of perturbation effectively confines the two 
hyperfine states in separate axisymmetric harmonic potentials with the same trap 
frequency uj_. The trap centers are offset spatially in the x-y plane by a distance 
r from the center and rotate about the symmetry axis at angular frequency 0'. 

The underlying physics for creating a vortex can be understood by considering 
the co-rotating frame with the trap centers at the angular frequency Q' of the 
perturbation H\. In this frame, the energy of a vortex with one unit of angular 
momentum is shifted by Ml' relative to its value in the laboratory frame. When 
this energy shift is compensated for by the sum of the detuning energy HA and the 
small chemical potential difference between vortex and non- vortex states, a resonant 
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Fig. 1. (a) The experimental technique used to create a vortex at JILA. An off-resonant laser 
provides a rotating gradient in the AC Stark shift across the condensate as a microwave that 
causes detuning is applied, (b) The level diagram shows the microwave transition to very near the 
1 2} state and the modulation due to the laser rotation frequency that couples only to the angular 
momentum i = 1 state when u> = A plus a small chemical potential difference between the vortex 
(£ = 1) and nonvortex (I = 0) state, where u> corresponds to f2' in the text. In the figure, the 
energy splitting (< 1 Hz) between the 1 = 1 and 1 = states is exaggerated. [Matthews et al, 
Phys. Rev. Lett. 83 2498 (1999), reproduced with permission. Copyright (1999) by the American 
Physical Society.] 

transfer of population can occur. From an initial nonrotating component |1), the |2) 
component is transferred to the state with unit angular momentum by controlling 
precisely the time of turning off of the coupling drive; this scheme was confirmed 
to work by numerical simulation of Eqs. (31) and (32) 12 . A typical structure shows 
that the |2) component has a vortex with a single quantized circulation at the 
center, whereas the nonrotating |1) component is located at the center and provides 
a "pinning" potential that stabilizes the vortex core [see Figs. 3 and 4]. 

Following the above method, Matthews et al. created a vortex for the first time 
in trapped BECs 11 . To realize the above configuration, they used a microwave and 
a movable laser beam with a spatially inhomogeneous profile on the condensate 11 . 
The laser beam was rotated around the initial condensate as in Fig. 1(a), which 
reproduces the perturbation given by Eq. (33). The resonant transfer of popula- 
tion from the nonrotating condensate into the vortex state was accomplished by 
detuning the microwave frequency and rotating the laser beam at an appropriate 
frequency for the resonant coupling. As shown in Fig. 1(b) for large detunings, the 
energy resonance condition means that atoms can change only the internal state 
through the coupling of the time-varying perturbation, and therefore must follow 
any selection rule that the spatial symmetry of that perturbation might impose. It 
should be emphasized that it is not simply a mechanical force of the optical field 
that excites the vortex. If one changes the sign of detuning A at a fixed trap rotation 
frequency, a vortex with opposite circulation will be created. Vortices with oppo- 



2, 2008 15:1 WSPC/INSTRUCTION FILE 2vorrev 



14 K. Kasamatsu, M Tsubota and M. Veda 

site circulations experience opposite energy shifts in transforming to the rotating 
frame, and therefore they must be detuned in an opposite sense to obtain resonant 
coupling. 

Experimentally it is possible to put the initial condensate into either the |1) or 
1 2) state, and then make a vortex in the |2) or |1) state. In the case of 87 Rb atoms, 
the intracomponent and intercomponent scattering lengths are in the proportion 
ai : ai2 : a,2 = 1.03 : 1 : 0.97 with the average of the three being 55 A. For the case 
without vorticity, the atoms with the larger scattering length a\ in the state |1) 
formed a lower-density shell around the atoms with the smaller scattering length 
a,2 70 [see Fig. 3(a)]. The stability properties of the vortex states are strongly depen- 
dent on the relation of these scattering lengths. The time evolution of the vortex 
can be observed over time scales from milliseconds to seconds in the experiment. In 
Ref. 11, the vortex was stable in only the state with the vortex in the |1) state with 
the |2) state in the core [Fig. 2(a,b)]. The other state [Fig. 2(c)], in which the vortex 
was in the |2) state, showed an instability with the |2) vortex sinking in toward the 
trap center and breaking up. Further discussion of dynamical stability is in Sec. 4.5. 



4.2. Structure of axisymmetric vortex states 
4.2.1. Thomas-Fermi approximation 

First, we consider the equilibrium structures of single vortex states in two- 
component BECs. As observed in Ref. 11, the simplest vortex structure consists 
of one circulating component that surrounds the other nonrotating component. 
Here, both wave functions have axisymmetric profiles. The authors in Ref. 78, 80, 
81 used the Thomas- Fermi approximation to extensively analyze this system. Here 
we summarize their findings. 

The starting point for the calculation is the coupled GP equations (21) and (22). 
Let us assume that the condensate wave functions have one quantized vortex at the 
center as \&i(r) = fi(r, z)e %liB ; then Eqs. (21) and (22) become 



& 



2 



2mi \ dr 2 

h 2 ( d 2 



2to2 \ dr 2 



1 d 


d 2 


_ V 


r dr 


dz 2 




1 d 


d 2 


72 


r dr 


*" dz 2 


r 2 



+ V 2 + 52/52 + .912/51 - 72^ 



/i = (34) 
h = M2/2, (35) 



where pi = ff is the condensate density and the trapping potential is given by 
V{ = mi(oj 2 ± r 2 + lo 2 z z 2 )/2. The axisymmetric vortex state is obtained by solving 
Eqs. (34) and (35) with (71,72) = (1,0) or (0, 1). 

Under the Thomas-Fermi approximation, the quantum pressure terms (deriva- 
tives with respective to the coordinates r and z) in Eqs. (34) and (35) are neglected. 
Then, the density profiles of the condensates can be calculated in a manner similar 
to that for the non-vortex case shown in Sec. 3. In the region where only one wave 
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Fig. 2. (a,b) Two examples of the free evolution of a magnetically trapped vortex in the |1) state. 
The lifetime of this configuration was much longer than the trap oscillation period (128 ms). (c) 
The free evolution of a vortex in the \2) state. This vortex exhibits unstable dynamics. The vortex 
first contracts and then splits into fragments. Each column is obtained from a single run of the 
experiment, where the origin of time t was chosen to be the time when the vortex was created (t 
is common to every row). The |1) and |2) state images appear different due to different signs of 
probe detuning. [Matthews et al, Phys. Rev. Lett. 83 2451 (1999), reproduced with permission. 
Copyright (1999) by the American Physical Society.] 



function is nonvanishing (pj ^ and pk = for j ^ k), Eqs. (34) and (35) are 
decoupled and the solution is 

V t - h 2 ll /{2m l r 2 ) 



Pi 



Hi 



(i = l,2). 



(36) 



On the other hand, in the region where the wave functions overlap (fi > 0), one 
obtains 



Pi 



P2 



1 



5i r i2 
1 

92^12 



Pi-Vi 



M2 - V 2 



2mi r 2 

fi 2 72 
Irrii r 2 



512 
.92 
912 

9i 



'■ — M2 - V 2 



— \Hi-Vi 



72 



2m,2 r 2 
h 2 



7i 



2r7Ji r 2 



(37) 



(38) 



Here, the chemical potential has been rewritten as pi — jiQ — > pi. The term 
h 2 -fi/(2mir 2 ) represents the density singularity of the vortex core. To obtain the 
density profile, we have to calculate the chemical potentials pi from the normaliza- 
tion condition j d 3 rpi = TVj. This leads to a problem involving two coupled integral 
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Fig. 3. Particle density of a vortex state in a trapping potential with parameters a\ = 5.67 nm, 
a2 = 5.34 nm, a\2 = 5.50 nm, uij_ = 7.8Hz, and ui z /w± = l/%/8, and Ni = N2 = 5 X 10 s , taken 
from the JILA experiments. The wave function is normalized as J d 3 r\^fi\ 2 = 1. The contour plots 
in (a)-(c) show the density profile pi = \"$!i\ 2 (i = 1,2) in the (r, z) plane obtained using the 
Thomas-Fermi approximation (top) and numerical simulations (middle) for (a) (7i,72)=(0,0), (b) 
(1,0) and (c) (0,1). The white (black) region represents the high (low) density region. The bottom 
plots show the corresponding radial density profiles p\ (solid curve) and p2 (dashed curve) for 
the equilibrium state at z = 0. The thin and bold curve represent the results obtained with the 
Thomas-Fermi approximation and numerical simulations, respectively. 

equations that can be solved numerically. Choosing the parameters corresponding 
to those of the JILA experiment 11 , namely, a\ — 5.67 nm, 02 = 5.34 nm, 012 = 5.50 
nm, uj± = 7.8Hz, uj z /uj± = 1/V&, and Ni = N 2 = 5 x 10 5 , we calculate the density 
profile for (71,72) = (0,0) as well as (1,0) and (0,1) by both using the Thomas- 
Fermi approximation and numerical simulation, and the results are shown in Fig. 3. 
When neither component has nonzero vorticity [Fig. 3(a)], the components phase 
separate due to the intercomponent repulsion. Specifically, to decrease the intra- 
componcnt mean field energy, the component with the larger intracomponcnt 
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scattering length forms a shell outside the 5*2 component in the central region. 
The density contours in the (r, z) plane show that |^i| 2 extends more along the 
z axis than the r axis. When the \&i component has a vortex, a density hole as- 
sociated with the vortex core appears in the density of and penetrates through 
the z-axis as shown in Fig. 3(b). Then, the ^ 2 component fills the vortex core 
because of the intercomponent repulsion, forming a "cigar-shaped" profile. As a re- 
sult, the characteristic scale of the core size becomes larger than the healing length 
£ = h/(2mgipi), which characterizes the core size of a single-component BEC. In 
this case, the centrifugal force associated with the vortex makes the component 
expand radially, which allows a decrease in the intracomponent mean-field energy 
of the \&i component rather than that of the nonvortex state. On the other hand, 
when the ^2 component with the smaller intracomponent scattering length has a 
vortex, the structure is different from the (71, 72) = (1,0) solution as shown in Fig. 
3(c). First, the core size is smaller than that of the (1,0) solution. Second, while a 
fraction of the 'J'i component fills the vortex core, the rotating '5 2 component is sur- 
rounded by a shell of the excessive ^1 component, despite the centrifugal force. As 
discussed below, the latter configuration is dynamically unstable 82 ' 83 ' 84 . The above 
consideration shows that the Thomas-Fermi approximation captures some quali- 
tative features of the solutions 78,80,81 . However, one cannot neglect the quantum 
pressure term for the description of the vortex states in this system; in particu- 
lar, the quantum pressure is necessary to determine the boundary profile between 
segregated components 136,137 . 



4.2.2. Pseudospin representation: Nonlinear sigma model 

To analyze the vortex structure, we present a 2-D (two-dimensional) model 92 that is 
based on the concept of "pseudospin" 79 > 91 > n o,i5o _ ^ ne S p mor order parameter of two- 
component BECs allows us to analyze this system as a spin-1/2 BEC 79,93,94,110,150 . 
To simplify the problem, we assume that all atoms have the same mass (mi = m2 — 
m) and the trapping potentials are equal (Vi = V2 = V). We map the Hilbert space 
of the system into pseudospin space by introducing the normalized complex-valued 
spinor x(r) = [xi( r )> X2(r)] T = [| Xi I e* 01 , |X2|e j6l2 ] T . Also, we decompose the wave 
function as — y PT(r)xi( r )- Here, pt = Pi + P2 is the total density, and the 
spinor thus satisfies 

lxi| 2 + lx 2 | 2 -l. (39) 

In terms of the spinor, the spin density is defined as S = x(r)cr x( r ); where <r is the 
Pauli matrix. Explicit expressions of S = (S x , S y , S z ) read 

S x = (X1X2 + X2X1) = 2|xi||X2| cos(0i - 2 ), (40) 
S y = -i( X *iX2 - X2X1) = -2|xi||X2| sin(0i - 62), (41) 
S z = (|xi| 2 -|X2| 2 ), (42) 
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where the modulus of the total spin is |S| = 1 everywhere. We assume that *i and 
\&2 represent the up and down components of the spin- 1/2 spinor, respectively. The 
nonzero spin projection on the x-y plane implies a relative phase coherence between 
the up and down components. 

The transformation of Eq. (14) to the pscudospin representation is 
straightforward 92 . The result is 

^ + ^(VSf + ^(v cff -Oxrf 

(43) 



E= d 



1?>, 



+Vp T + ^(c + c 1 S z + c 2 S 2 z ) 



where we define new coupling constants as 

9i + 92 + 2gi2 



co 



■1 



Cl = (44) 



c 2 = 



.91 + 92 - 2.9 



12 



and an effective velocity field as 



Veff = ^(Xi V Xi _ XiVxt + X*2^X2 - X2VX2) 



2^ 



yQ _|_ S_z(SyVS x S x VSy) 



(45) 



which depends on the gradient of the total phase 6 = Q\ + 9 2 and that of the 
pseudospin. 

The form of Eq. (43) is analogous to the classical nonlinear sigma model (NLctM) 
for Heisenberg ferromagnets in which only the (VS) 2 term appears. A NLuM al- 
lows one to study a various type of "topological excitations" such as monopoles, 
domain walls, and skyrmions 155 . Because we have changed merely the representa- 
tion of the Gross-Pitaevskii energy functional, Eq. (43) is still an exact description 
of the two-component BEC. The four degrees of freedom of the original condensate 
wave functions ^1 and \&2 (their amplitudes and phases) have now been expressed 
in terms of the total density px, the total phase 6, and two of the spin density 
components (S x , S y , S z ) (one of which is determined besides its sign by the other 
two because of |S| = 1). 

Unlike the classical NLcrM case, the system here has several unique features 
that are revealed in Eq. (43): (i) the total density pr, which is a prefactor of the 
(VS) 2 term and effectively 'stiffens' the pseudospin, is spatially inhomogeneous be- 
cause of the trapping potential; (ii) there is a hydrodynamic kinetic-energy term 
wpT(v c ff — H x r) 2 /2, which is associated with the topological excitations; (iii) 
if 9i 92 7^ .912, the anisotropic terms involving the coefficients c\ and c 2 effec- 
tively break the "SU(2)-invariance" 79 ' 92 ' 93 ' 156 because then the rotation of S re- 
moves the energy degeneracy. The coefficient c\ can be interpreted as a longitudinal 
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(pseudo)magnetic field that tends to align the spin along the z-axis. The term with 
the coefficient c 2 may determine the spin-spin interaction associated with S z , which 
is antifcrromagnctic for c 2 > and ferromagnetic for c 2 < 88 . If mi ^ m 2 or 
Vi j^Vi, other anisotropic terms appear in Eq. (43). 

In the following, we confine ourselves to the 2-D problem for simplicity. Two- 
dimensional calculations of the GP equations (21) and (22) are carried out by sepa- 
rating the degrees of freedom of the original wave function as ^(r) = tpi( x i y)4>( z )- 
This separation yields the dimensionless 2-D coupling constants m — 4ira,ir]N and 
U12 = 47rai2?7-/V with t) = J dz\(j>(z) | 4 / / dz\(j)(z)\ 2 [see Ref.37 for more details], 
where N — N\ + N 2 is the total particle number. The dimensionless GP equations 
become 



1 



~2 V +y +ui\ipi\ 2 + u 12 \ip 2 \-nL z )Vi =MiVi, 
1 r 2 

-;:V 2 + — + U 2 \lp 2 \ 2 + Ul2|^l| 2 - QL X ) i> 2 = (12^2 



(46) 
(47) 



where the length and time are measured in units of &ho = \J h/muj± and w^ 1 , which 
are the characteristic scales of the radial trapping potential V — moj 2 _r 2 /2. Also, 
the NLaM (43) becomes 



E 



^(Vv^) 2 + ?(VS) 2 + ^(v cff - SI x r) 2 



Pt, 



-yPT + ^ (co + c x S z + c 2 S z ) 



(48) 



with c = (ui + u 2 + 2mi 2 )/4, c\ = (u\ - u 2 )/2, and c 2 = (ui + u 2 - 2m 2 )/4. 



4.2.3. A skyrmion: Numerical analysis 

First, we consider the axisymmetric vortex state for the case gi = g 2 = g\ 2 (c.\ = 
c 2 = 0), in which case the system possesses SU(2) symmetry. In the pseudospin 
picture, the axisymmetric vortex observed in Rcf. 11 can be interpreted as a type 
of spin texture called a " skyrmion" 79 ' lw . Typical examples are shown in Fig. 4. In 
Fig. 4(a), we assume that the -01 component has one singly-quantized vortex at the 
center of the trap. At the center of the atomic cloud, the ipi component vanishes, 
so that the pseudospin points down according to the definition of the spin S z in Eq. 
(42). The spin aligns with a hyperbolic distribution as {S x ,S y ) oc (— x,y) around 
the singularity at the center. At the edge of the atomic cloud, the ?A 2 component 
vanishes, and the pseudospin points up. In between, the pseudospin continuously 
rolls from down to up as shown in the bottom of Fig. 4. This "cross-disgyration" 
texture is often referred to as a skyrmion in analogy to the work of Skyrme 6 . In the 
field of superfluid 3 He, this skyrmion texture is referred to as an Anderson- Toulouse 
vortex 79 . If ip 2 component has a vortex, the texture exhibits a "radial-disgyration" as 
illustrated in the middle panel of Fig. 4(b), in which the spin in the x-y plane aligns 
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4 







Fig. 4. Results of simulations of 2-D coreless vortex states. Top panels: the density profiles con- 
sisting of (a) the rotating tjix component and the nonrotating ip2 component and (b) the non- 
rotating ipi component and the rotating ip2 component. The parameters are ui = U2 = 1112 
(co = 1000, ci = C2 = 0) and Q = 0.15. For the calculation, we fixed only the total particle 
number N = Ni + N 2 . Then, the solution converges to N1/N2 = 2.465 in (a) and 0.406 in (b). 
Middle panels: vector plots of the spin texture projected onto the x-y plane corresponding to (a) 
(left) and (b) (right). The region is [—5 < x,y < 5]. The bottom panels: the cross-sections of the 
spin texture along the a;-axis at y = corresponding to the vector plots of the middle panels. 

as (S x ,S y ) oc (x,y) around the singularity. These two skyrmions are degenerate in 
the SU(2) symmetric case because they are interconverted by the overall spin; in 
fact, all configurations obtained by any global spin rotation are degenerate. 

The size of the skyrmion is determined by the interaction coefficients and the ro- 
tation frequency 92 , and is related to the ratio N\/N%. When we minimize the energy 
to determine the ground state, we fix the total particle number, not each particle 
number. This procedure gives us the true thermodynamically stable configuration of 
the skyrmion. Although each particle number is conserved in experiments performed 
to date (this case will be discussed below) , the population imbalance determined as 
above can be realized by controlling the strength and time of the coupling drive in 
the phase imprinting method. 

It is known that the skyrmion in a 2-D system has a topological invariant defined 

as 

Q=^- [ d 2 re ij S ■ <3;S x djS, (49) 
on J 
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Fig. 5. (a) Topological charge density q(r) defined by Eq. (50). (b) The effective velocity field v c fj 
defined by Eq. (45), corresponding to the solution shown in Fig. 4. The insets are cross-sections 
of q(r) and |v e ff | along the y = line within the range —5 < x < 5. 



which is known as the topological charge or the Pontryagian index 155 . The skyrmion 
with any spin profile is shown to have Q = ±1, whose sign depends on the direction 
of the circulation. The integrand of Eq. (49) is the topological charge density g(r) 
associated with the vorticity. The topological charge density can be derived from 
the effective velocity v g as 92 ' 110,111 : 

q(r) = i- e «S • OS x djS = ^(Vx v cff ) z , (50) 

where we used the relation J2iXi^Xi = ~SiXi^Xi (* = 1)2). The topological 
charge density q(r) characterizes the spatial distribution of the skyrmion. Figure 5 
shows the spatial distribution of q(r) and the corresponding v e fi -field for the solution 
of Fig. 4(a). The topological charge is distributed around the center and, contrary to 
the case of a conventional vortex in a single-component condensate, |v e ff| vanishes 
at the center. This makes a coreless vortex (skyrmion) without a density dip in the 
total density. 



4.2.4. A skyrmion: Variational analysis 

To study the physical properties of the skyrmion in more detail, we used a variational 
analysis on the NLctM described by Eq. (48). In the following, we confine ourselves 
to the case in which the ipi component has a vortex. The original NLctM, where 
only the (VS) 2 term is present in Eq. (48), admits a skyrmion solution whose 
explicit analytic expression is known 155 . Here we seek for a more general form of 
the skyrmion solution of Eq. (48) as a variational function, which can reproduce 
well the numerical solution. The skyrmion solution of Fig. 4(a) may be parametrized 
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as 157 



AXxe- ar2 / 2 ~AXye- ar2 / 2 _ r 2 - 4A 2 e -° r2 

S * ~ r 2 + 4A 2 e -ar- A f ~ r 2 + 4A 2 e -ar" 6 * ~ r 2 + 4A 2 e -aH ^ 

with |S| = 1. Equations (51) with a = correspond to the explicit skyrmion 
solution of the classical NLctM 155 . The variational parameters A and a determine 
the size and the shape of the skyrmion. Typically A represents a characteristic 
size of the domain over which the direction of the spin is maintained and a is the 
degree of spatial variation of the spin inversion. If we take both A and a as variable 
parameters, the number of particles in each pseudospin component is not conserved, 
but the total particle number is conserved. In the limit A — > (oo) with fixed a all 
spins point up (down) along the z-axis, which means a perfect polarization of the 
particle number such as N —> N\ (A^). Alternately, for fixed A, the limit N — > N\ 
(N2) corresponds to a — > +00 (—00). Because the ipi component has one singly- 
quantized vortex at the center, the total phase 9 in Eq. (48) is 

9 = tan" 1 -. (52) 



Substituting Eqs. (51) and (52) into Eq. (48) yields 



E = j d 2 rji(V v /^) 2 + Wt + ^ 



co + ci 



4AV 



r 2 + 4A 2 e- Qr2 



+ C 2 



+ 4A 2 e 



(53) 



Thus, there are now fewer degrees of freedom in the energy functional: the total 
density px and the two variational parameters A and a. In the above equation, we 
introduced an effective, radially symmetric, confining potential 



r 2 + A\ 2 e- ar ' 2 [{ar 2 + l) 2 + l] fir 2 



r 



2 



VcS ~ 2(r 2 + 4A 2 e-«'' 2 ) 2 r 2 + 4A 2 e -«'' 2 + 2 ' (54) 

which determines the shape of the total density px- The total density px can be 
calculated by solving the equation 



(v 2 Vp?) 

2v^ 



V'eff + PT 



c o + ci „ , _ —5 + c 2 



r 2 + 4A 2 e- Qr2 V?' 2 +4A 2 e 



2t 



= /i, (55) 



where the chemical potential \i is fixed by the normalization condition J rpx = 1- 
Using the condition J cPrpT = 1, we numerically calculate the chemical potential 
fi in Eq. (55) with the Thomas-Fermi approximation, where px is assumed to be 
zero in the region in which px is negative. Then, the total energy E becomes a 
function of (A, a) by using the Thomas-Fermi expression for px- When the size 
of the skyrmion is large enough 92 , the Thomas- Fermi approximation can be used 
to evaluate px, and then the energy can be minimized with respect to only two 
variational parameters A and a. Otherwise a density singularity associated with the 
vortex core appears in px so that the quantum pressure term cannot be ignored. 
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Fig. 6. (a) Optimized values of the variational parameters A (solid curve) and a (dashed curve) 
as a fuction of rotation frequency Q and for co = 1000, 2500, 10000 with ci = C2 = 0. (b) The 
cross-sections of \ipi\ 2 (solid curve), |i/"2| 2 (dashed curve) and the total density px (dotted curve) 
along the ai-axis at y = 0, where bold and thin curves represent the results obtained from the 
numerical and variational calculations, respectively. 



For the SU(2) symmetric case without rotation (fi = 0), the minimum energy- 
occurs when A — > oo and a — > — co, which implies that N = N2, i.e., complete 
polarization of the particle number. It follows that the global minimum for Q = 
is a nonvortex state for any value of N2 for a fixed total particle number N. This is 
because the interaction energy terms now satisfy the SU(2) symmetry so that the 
energy of the nonvortex state is degenerate under the change of the ratio N1/N2 
with fixed N. However, above a certain critical value of f2, there appears an energy 
minimum for particular finite values of A and a. Since the minimized energy is lower 
than that of the nonvortex state with A — > +00, an appropriate external rotation 
can provide global stability for the skyrmion. Figure 6(a) shows the values of A and a 
that give the minimum of the total energy for cq = 1000, 2500, and 10000. The size 
of the skyrmion decreases with increasing f2, as revealed by a decrease in A m i n . The 
same effect occurs for an increase in a. The energy minimum appears at a certain 
critical frequency of £1 where A m i n (a m i n ) diverges to +00 (—00) and the critical 
frequency decreases as cq increases. Therefore, the condensates with larger Co can 
have a stable skyrmion at a lower rotation frequency. The comparison with the 
numerical solution in Fig. 6(b) shows that the optimized variational functions given 
in Eqs. (51) and (52) nearly reproduce the exact numerical solution. Hence, our 
approach here greatly improves upon the analytic treatment for the vortex states 
beyond the usual Thomas-Fermi approximation in Refs. 78, 80, 81. However, full 
numerical calculations show that for cq = 1000, 2500, 10000 additional vortices are 
nucleated above f2 > 0.17, 0.11, 0.05 and the ansatz (51) for a single skyrmion is 
no longer valid. 
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The parameters c\ and c 2 in Eq. (53), the values of which can be controlled 
through the scattering lengths, predominantly determine the size of the skyrmion. 
The coefficient ci, which is proportional to the difference between the intracom- 
ponent interactions u\ and tt 2 , may be regarded as a (pseudo)magnetic field that 
tends to cause the spin to align along the z-axis; the positive (negative) c\ arranges 
the spins downward (upward). For the negative ci (ui < u 2 ) the stable size of the 
skyrmion shrinks as |ci| increases, which implies an increase in the fraction of the 
rotating ipi component and a decrease in the fraction of the ip 2 component at the 
vortex core, (i.e., the spins align upward). For the positive c\ {u\ > u 2 ) the value 
of A increases rapidly and eventually goes to infinity (concurrently, a goes to — oo), 
corresponding to the perfect "spin-down" alignment (occupation of the nonrotating 
ip 2 component only). 

The sign of c 2 determines the nature of the spin-spin interaction associated with 
S z . For positive c 2 (u\ + u 2 > 2ui 2 ) the interaction is antiferromagnetic, that is, a 
spatial mixture of the spin- up and spin-down components. The presence of a vortex 
causes an effective phase separation. Then, there appears to be no significant change 
of the spin profile compared to that of c 2 = in Fig. 4(a). In contrast, for negative 
c 2 {ui + u 2 < 2ui 2 ), the system enters the ferromagnetic phase, where the spin 
domains are spontaneously formed. When there is a skyrmion in the ipi component, 
its size shrinks with |c 2 | in a way similar to the |ci|-depcndcnce. This is because 
each particle number is not conserved in the present case; the energetically favorable 
configuration in a ferromagnetic phase tends to an overall spin polarization, that 
is, a perfect polarization of the particle number. In this regime, depending on the 
rotation frequency, there are two energy minima, one of which leads to a perfect 
polarization of the tpi vortex state, but the other leads to the ip 2 nonvortex state, 
which is stable. 

So far, we have allowed a change in the particle number of each component 
to facilitate the calculation of the minimum skyrmion size. However, the actual 
experiments on two-component BECs have been done under a condition in which 
each particle number is fixed. This restriction can be taken into account by imposing 
the relation 

m Jd^rp T (l + S z ) 

N 2 jd?rp T {l-S z y [00) 

For a given A the value of a is uniquely determined by Eq. (56), where both pt and 
S z are functions of two variational parameters. Therefore, the energy minimization 
can be done with respect to one variational parameter, which we choose to be A. Let 
us here consider the case N\/N 2 = 1, and investigate the stable size of a skyrmion. 
In this case, we find that the stable size is not affected by the change in rotation 
frequency Q. Figure 7 shows the optimized values of A and a as a function of c\ and 
c 2 . No perfect polarization of the particle number occurs, so that the skyrmion may 
exist for all values of c\ and c 2 ; moreover, in this case, even if the optimal variational 
parameter exist, they do not ensure a local minimum of the energy. In Fig. 7(a), for 
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(a) 3 



(b) 3 





c i=- 50 








Fig. 7. Optimized values of the variational parameters A (solid curve) and a (dashed curve) as 
a function of (a) c\ and (b) C2 under the condition of fixed particle number of each component 
Ni = N 2 for c = 10 3 , n = 0.15. We fix c 2 = for (a) and ci = for (b). The right panels show 
the cross sections of the optimized variational functions \ipi\ 2 (solid curve), \ip2\ 2 (dashed curve), 
and the total density px (dotted curve) along the y = line ((a) top: ci = —50, bottom: ci = 50 
and (b) top: C2 = —30, bottom: C2 = 30). 



positive ci, to enlarge the "spin-down" domain containing the ip2 component, one 
should increase both Amm an d c*min- When this domain increases, so does the size of 
the vortex core of the ipi component. For negative c±, the domain with the rotating 
ipi component tends to increase. As a result, the vortex core of the ipi component 
shrinks and the excess ^2 component, which fills the vortex core, protrudes into the 
outside of the ipi component. In the ci dependence, although there are no drastic 
change in the antifcrromagnetic ci > regime, for ci < a rapid increase occurs in 
A m i n and a m i n when \c-i \ increases [Fig. 7(b)]. This means that spin-up or spin-down 
domains grow and their boundary becomes sharper for larger |ca | • These results are 
consistent with the numerical ones shown in Fig. 3. 



2, 2008 15:1 WSPC/INSTRUCTION FILE 2vorrev 



26 K. Kasamatsu, M Tsubota and M. Veda 

4.3. Structure of nonaxisymmetric vortex states: a pair of merons 
4.3.1. A meron pair: Numerical analysis 

In this section, we discuss a nonaxisymmetric vortex state. The spin texture of 
this state involves a pair of "merons" 4 or "Mermin-Ho vortices" 79,98,110 ("meron 
(/iepocr)" means "fraction" in Latin). Our previous study revealed that a meron 
pair (a vortex molecule) can be stabilized thcrmodynamically in a two-component 
BEC if it is rotated and has an internal coherent coupling 91 . The two components 
interact not only through their mean-field interactions but also through the relative 
phase of the order parameters. Combination of these two effects enables us to explore 
a new regime of rich vortex structures which are absent in the conventional binary 
system 78,80,81,86,87,88 . If the strength of the coupling drive is increased gradually 
from zero and its frequency is adiabatically tuned to be resonant, one can obtain a 
stationary state with a nearly equal- weight superposition of the two states 150 . The 
internal coherent coupling can be included in the formulation through Eq. (30), and 
we obtain the time- independent, coupled GP equations by substituting &i(r, t) = 
<I/i(r)e -vt / fi into Eqs. (31) and (32). Both components are characterized by the 
common chemical potential /i because, in the presence of the coherent coupling, the 
conserved quantity is the total particle number rather than the individual particle 
numbers. For simplicity, the detuning parameter A is assumed to be zero. In terms of 
the pseudospin picture, the Rabi term given in Eq. (30) is represented as — u>hPtS x , 
which plays a role of a "transverse magnetic field" along the x-axis. 

A typical solution of the vortex state with coherent coupling is shown in Fig. 8. 
Each component has one off-axis vortex and the density peak of one component is 
located at the vortex core of the other component. This results in a coreless vortex 
in which the total density has no singularity. Because of the coherent coupling, the 
profile of the relative phase </>(r) = 61 — 62 plays an important role in optimizing the 
structure, which is shown in Fig. 8(b). The relative phase shows that the central 
region is characterized as a vortex-antivortex pair. In other words, the two vortices 
are connected by a branch cut of the relative phase with the 2tt difference, which 
is a characteristic domain wall structure in two-component BECs with internal 
coherent coupling 158 . Then, the vortex in one component attracts that in the other 
component, due to their being bound by the domain wall, thus forming a "vortex 
molecule" 91 . As ojr increases, the size of the pair decreases as seen in Fig. 9. Beyond 
cjr ~ 3.0, the two vortices completely overlap despite the intercomponent repulsive 
interaction. 

The corresponding spin texture is shown in Fig. 8(c) and (d). Because of the 
transverse "magnetic field" , the spins are oriented in the x direction everywhere ex- 
cept in the central domain- wall region where they tumble rapidly by 27r. There exist 
two points, corresponding to the vortex centers, at which S is parallel to the z-axis. 
In Fig. 8(c), the spin around the singularity with S = +z has a radial distribution, 
characterized by (S X7 S y ) oc (— x, —y), and rolls by 90° as it goes outward, becoming 
perpendicular to the x-axis. For the spin around S = — z, the distribution is instead 
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y 

Fig. 8. A simulated meron pair. The values of the parameters are ui = U2 = u\2 = 1000 
(co = 1000, c\ = C2 = 0), f2 = 0.15 and u;r = 0.05. (a) Density profile of the condensates l^ll 2 
and \il>2\ 2 for Co = 10 3 , c\ = C2 = 0, Q = 0.15 and ojpi = 0.05. The two components have the same 
particle number, (b) Gray scale plot of the relative phase if) = 9± — 02 = arg?/>i — arg?/>2- Arrows 
show the direction of the circulation in relative-phase space, (c) Vectorial representation of the 
pseudospin texture of the vortex state (a), projected onto the x-y plane, (d) Cross section of the 
spin texture (c) along the x axis at y = 0. 

hyperbolic (S Xl S y ) cx (x,—y). This texture is known as a "radial-hyperbolic" pair 
of merons, which has been discussed in the study of topological defects in superfluid 
3 He (Ref. 2) and a double-layer quantum Hall system 4 . Since the trapping poten- 
tial is axisymmctric, the energy is degenerate with respect to the orientation of the 
molecular polarization. If the molecule is polarized along the y-axis, the texture 
features a "circular-hyperbolic" pair, but it can continuously transform into the 
radial-hyperbolic pair. 

We find that, in the SU(2) symmetric case, i.e., u\ = U2 — U12 (ci = C2 = 0), 
both the topological charge density q(r) and the effective velocity v c ff exhibit a 
radially isotropic profiles like that in Fig. 5, even though each component forms a 
nonaxisymmetric vortex configuration. This is due to the fact that the meron pair 
and the axisymmetric skyrmion are equivalent as a topological excitation for the 
case of ui = U2 = U12; they can transform to each other via an overall rotation of 
the pseudospin. Then, giving an infinitesimal value to o;r (i.e. applying a transverse 
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Fig. 9. Cross sections of a meron pair (|i/>i| 2 : solid-curve, IV'al 2 : dashed-curve, and px'- dotted- 
curve) along the x-axis at y = for cq = 10 3 , ci = C2 = 0, Q = 0.15 and =0, 0.02, 0.2, 0.8, 2, 
3. 



magnetic field) is sufficient to stabilize a nonaxisymmctric configuration. We thus 
obtain a further insight of the nonaxisymmetric vortex by rotating the basis of the 
spinor so that the internal coupling becomes simpler. The Rabi term plays the role 
of a transverse "magnetic field" , making the x-axis as a preferred axis. Actually, by 
rotating the spinors as ip± = (ipi ± i/'2)/v / 2 (the basis along with pseudo-cc axis), 
the coupled GP equations become 



1 /V 



2 V i 



fl x r + V + CqPt 



2 V i 



V + cqPt 



i/j- = (p- Wr)-0— 



(57) 
(58) 



Here, the internal coupling is just the chemical potential difference between the "+" 
and "— " components. Then, the nonaxisymmetric structure in Fig. 8 is transformed 
to the axisymmetric vortex state in Fig. 4, that is, a skyrmion, where, the vortex 
core of the "+" component is filled with the nonrotating "— " component. As one 
increases ur the number of the "— " particle decreases due to a decrease in the 
chemical potential of the ip- component, and thus the vortex cores eventually be- 
come empty. This corresponds to the overlap of the vortex cores as shown in Fig. 9. 
Therefore, the Rabi term alone does not break the axisymmetry of the topological 
excitation. However, inclusion of both the Rabi term and unequal coupling con- 
stants u\ ^ U2 ^ 1*12 induces an anisotropy of the meron pair as will be discussed 
later. 



I 



I 
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4.3.2. A meron pair: Variational analysis 



In the case of u\ = u 2 = u\ 2 {c\ = c 2 = 0), the spin profile of the meron pair in 
Fig. 8 may be parametrized as 91 ' 92 

_ r 2 - 4X 2 e- ar2 -4Xye- ar2 / 2 -AXxe^ 2 ' 2 

b *~ r 2 + 4A 2 e -ar- b V ~ r 2 + 4A 2 e -ar" b * ~ r 2 + ^-ar* " ^ 9 > 

Here, we assume that the meron pair is polarized along the x-axis. The only differ- 
ence between this case and that for the skyrmion in Eq. (51) is that the forms of S x 
and S z have been exchanged. In this case, the ratio of the particle number given in 
Eq. (56) is always unity, because S z is an odd function. The locations of the vortex 
cores are determined by two extremes of S z , and are given by 

x 2 =4\ 2 e- ax \ (60) 

Then, it is natural to use the following form for the total phase 0: 

9 = tan " a -2Ae-^/2 +tan " 1 x + 2Xe~^ - (61) 

Substituting Eq. (59) and Eq. (61) into Eq. (48) with the Rabi term — ujrPtS X} wc 
obtain the following total energy: 



E 



-I*' 



I r — 4A e~ ar p 



(62) 



where the effective confining potential V e & is the same as that in Eq. (54). Because 
for any values of A and a Eq. (62) gives the same energy as that of the axisymmet- 
ric skyrmion given by Eq. (53), the skyrmion and the meron pair have the same 
optimized values of A and a if u>r = (i.e., for the SU(2) symmetric case), and their 
energies are degenerate. 

Just as the Ci-term may be regarded as a magnetic field along the z-axis, the 
Rabi frequency wr may be regarded as a magnetic field along the x-axis. Thus, the 
magnitude of cjr also influcences the stable size of the meron pair. The difference 
from the Ci-term is that the Rabi term is proportional to px instead of p\. We cal- 
culated the minimum values of A and a as a function of wr under a slow rotation. 
The result is shown in Fig. 10(a). As cjr increases, the minimized values of A (a) 
becomes smaller (larger), and eventually vanishes (diverges). This behavior corre- 
sponds to a decrease in the size d m of the meron pair (see in Fig. 9 for the definition) 
as shown in Fig. 10(b), where d m is calculated from Eq. (60). This indicates that 
the binding of the meron pair becomes stronger with increasing wr. The variational 
result agrees well with that obtained with direct numerical simulations. 

The meron pair arises from competion between repulsive and attractive 
interactions 91 . For well-separated merons, the repulsive interaction between them 
originates from the second and third terms of Eq. (48), which are the gradient en- 
ergy of the pseudospin and the hydrodynamic kinetic energy of the v e ff -field. The 
contribution of the other terms are nearly independent of meron separation except 
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Fig. 10. Properties of a meron pair, (a) Optimized values of the variational parameters A (solid 
curve) and a (dashed curve) as a function of for Co = 10 3 , CI = 0.15, c\ = 0, and C2 = 0. (b) The 
corresponding size 2d m of the meron pair (see Fig. 9 for definition) as a function of Wr, calculated 
from Eq. (60) with cq = 10 3 . Solid circles are the result obtained by numerical simulations. 



for small d m (< 0.30 for Co = 10 3 ). Then, the Thomas-Fermi approximation can- 
not apply to the evaluation of the total energy; A m ; n drops suddenly to zero with 
increasing cjr because there is no energy barrier associated with the quantum pres- 
sure of pt- On the other hand, the attractive force between the merons arises from 
a tension T<j in the domain wall of the relative phase, which is estimated to be 
~ Tdd m ; for a homogeneous system Td = 8|V'i| 2 |V'2| 2 A/pt with the characteristic 
domain size fc _1 = (|V'i||^'2|/2wr/9t) 1 ^ 2 (Ref. 158), and thus T<j oc s/lJb,. Thus, the 
competition between the repulsive force and the attractive force creates an energy 
minimum and the two vortices can form a bound pair. 



4.3.3. Effect of symmetry breaking terms 

Let us consider the effect of the symmetry breaking terms on a skyrmion and a 
meron pair. When m ^ U2 7^ U12 (ci, ci ^ 0), the axisymmetry of the topological 
excitation is broken. First, we neglect the spin-spin interaction term (ci term) and 
investigate the dependence of the stable structure on c\ and wr. This is equivalent to 
applying a longitudinal and transverse (pseudo)magnetic field. We prepare a stable 
axisymmetric skyrmion by applying the longitudinal magnetic field c\ = —10, and 
then turn on the transverse magnetic field wr. Figures ll(a)-(d) show variations 
of the spin texture as well as those of the cross-section of the condensate density 
along the y — line when the Rabi frequency wr is increased. The transverse 
magnetic field shifts the skyrmion off-center and converts it into a meron. At the 
same time, another meron enters from outside and forms an "asymmetric" meron 
pair shown in Figs. 11(c) and (d). This is a second-order phase transition because 
there is no energy barrier to destabilize the axisymmetric skyrmion with respect to 
the transverse magnetic field. In this case, the distribution of the topological charge 
q(r) moves from the center as increases, and there is no dramatic change in the 
global shape from the isotropic skyrmion of Fig. 5. As wr increases further, the 
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Fig. 11. Equilibrium structures of the spin texture in the presence of a longitudinal magnetic field 
(ci) and a transverse magnetic field (ur). The range of the plot is [—4 < x,y < +4]. The values 
of the parameters are f2 = 0.15, cq = 1000, c\ = —10, C2 = 0, and (a) ojr = 0, (b) ajR = 0.02, 
(c) ujr = 0.05, (d) ur = 0.10. The insets show the cross sections of IV'll 2 (solid curve), |V'2| 2 
(dashed curve), and the total density px (dotted curve) along the y = line within the range 
[-8 <x< +8]. 

peak of q(r) moves back to the center and its value goes to infinity when the two 
merons merge. 

Next, we study the effects of the spin-spin interaction term (c2 term) and the 
dependence of the stable structure on cjr (ci is fixed at zero). First, we show that, 
in contrast to the case of c\ , the structure of the meron pair is sensitive to C2 . Figure 
12 shows the equilibrium structure of the condensate density, the spin texture and 
the topological charge density q(r) for both the antiferromagnetic case [pz = —20) 
and the ferromagnetic case (02 = 20). For the antiferromagnetic case there is no 
significant difference in the density and spin profile, compared with the solution of 
C2 = in Figs. 8 and 9. However, the topological charge density (i.e., vorticity) 
is distributed anisotropically in such a way that its distribution is elongated along 
the direction of polarization of the meron pair. For the ferromagnetic case, spin 
domains are formed, which gives rise to a considerable change in the density and 
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Fig. 12. Equilibrium structures of the spin texture in the presence of a spin-spin interaction C2 
and a transverse magnetic field (iiij). The values of the parameters are fi = 0.15, Co = 1000, c\ = 0, 
lur = 0.05, and (a) C2 = 20, (b) C2 = —20. (top) The spin textures are displayed in the region 
[—4 < x, y < 4]. The insets show the cross sections of (solid curve), \ip2\ 2 (dashed curve), 

and the total density px (dotted curve) along the y = line within the range [—8 < x < +8]. The 
surface plots at the bottom are the spatial distributions of the topological charge density q(r). 

spin profile as seen in Fig. 12(b). Most spins align up or down along the z-axis 
on either side, sandwiching the domain wall across which the spin flips rapidly. 
If uj-h is increased further, the spins on both sides lay along the x-axis, forming 
the well-defined meron pair as in Fig. 9. Then, the topological charge density is 
distributed with the direction perpendicular to the direction of polarization of the 
meron pair, being concentrated on the domain-wall region. This anisotropy of the 
meron pair gives an interesting situation when the condensates undergo a rapid 
rotation; the anisotropic interaction between the meron pairs generates a distorted 
lattice of "vortex molecules" 91 , which will be discussed in Sec. 5.3. 

4.4. Topologically nontrivial skyrmions in a 3-D configuration 

In contrast to the topologically trivial skyrmion discussed above, we review here the 
topologically nontrivial skyrmion excitations in multicomponcnt BECs 93 ~ 96 . These 
skyrmions have complicated 3-D structures and are characterized by spin inversion 
in a finite region of space and a typical topological charge Q = 1. It is known that 
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topological considerations allow for these excitations, but it is not known whether 
or not such a configuration is energetically favorable. 

The first example of a nontrivial skyrmion was discussed by Al Khawaja and 
Stoof 93 . They considered uniform two-component condensates subject to the SU(2) 
symmetry and the radially symmetric profile of the skyrmion by assuming that the 
normalized spinor is characterized by x(r) = cxp{— 2iO(r) • S}\z, where S = er/2 
with the Pauli matrix er, Xz = (1,0) T , and fi(r) = u(r)r/r which determines 
a specific spin texture of the skyrmion. The physical meaning of this formula is 
that the average spin at a position r is rotated by an angle 2f2(r) from its initial 
orientation with the axis of rotation fi(r)/Q(r). The boundary conditions for oj(r) 
that represents the spin distribution for the skyrmion are the following. First, at 
r — > oo all spins must be oriented as in the ground state, that is, u{r — > oo) = 0. 
Along the z-axis the spins are also not rotated due to our assumption about f2; thus, 
to have a nonsingular texture of the spinor with a non-zero winding number, we 
must require w(0) = 2ir. Finally, to avoid the singular behavior of fi(r) itself, we use 
only functions u{r) with zero slope at the origin. Therefore, u)(r) is a monotonically 
decreasing function that starts from 27r at the origin and reaches zero as r — > oo. 

Al Khawaja and Stoof 93 used a variational ansatz ui(r) = 4cot~ 1 (r/A) 2 that 
satisfies the above boundary conditions, where A represents the size of the skyrmion. 
Substituting this ansatz into the energy functional in Eq. (14) with Vi = and £1 = 
gives the total energy as a function of the total density pt and A. For the given 
aj(r), they calculated pr{r) exactly by solving numerically the differential equation 
for pr(r) obtained by varying E[pr(r),X] with respect to pr{r). By substituting 
this density profile back into the energy functional, the energy of the skyrmion is 
expressed as a function of only A and the equilibrium properties can be obtained 
by minimizing this energy with respect to A. The resulting spin texture shows a 
complicated spin distribution; the spin (S z ) flips five times as it goes outward from 
the center in the x-y plane, and (S x ) and (S y ) exhibit quadrupole profiles with one 
radial node in the y-z plane 93 . 

A crucial problem is the stability of the skyrmion. In equilibrium, the energy is 
minimized for the vanishing size of the skyrmion 93 . However, for sufficiently small 
skyrmion sizes, a nonequilibrium stability mechanism starts to work. A number of 
atoms in the center of the skyrmion, hereafter referred to as core atoms, will be 
trapped by an effective 3-D potential barrier, that is, a repulsive shell with a finite 
radius that is created by the gradient term of the spin texture ft 2 |Vx(r)| 2 /2m. As 
the skyrmion shrinks in size, the barrier height of the repulsive shell increases and 
the radius decreases. This leads to a squeezing of the core, thereby increasing its 
energy, and ultimately stabilizing the skyrmion. In this case, it is not an equilibrium 
state of the condensate because the core atoms will tunnel through the barrier and 
give the skyrmion a finite lifetime. Detailed calculations of this tunneling rate 93 ' 159 
showed that only a few atoms in the core of the skyrmion arc needed to stabilize it 
and to give it a sufficiently long lifetime. 

Other interesting structures of nontrivial 3-D skyrmion were investigated in Refs. 
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94, 95, 96. The authors in Refs. 94, 95, 96 also considered the rotational operator 
of the spinor as x( r ) = exp{iA(r)er • (r/r)}x z , where A is a monotonic function 
imposed by the boundary condition for A as A(0) = and A = ir at the boundary of 
the condensate; for this boundary condition the topological charge becomes Q = 1. 
The corresponding condensate wave functions are 

(ipi{r)\ r—( cos[A(r)] - i sin[A(r)] cos 9 \ . . 

\M T )J \ -isin[A(r)]sin0exp(i0)y ' { ' 

Here, (A, 9, 4>) can be understood as the spherical angles of the 3-sphere (in four- 
dimensional parameter space) . Because of the topological stability of the skyrmion, 
any continuous deformation of Eq. (63), without altering the asymptotic boundary 
values, is still a skyrmion with Q = 1. With this realization of the skyrmion, \4>2\ 2 
vanishes both on the z-axis and at infinity, and is concentrated in a toroidal region. 
On the other hand, iV'il 2 vanishes on the circle 9 = ir/2 and r — X~ 1 (tt/2), a 
region where 1^2 1 2 is concentrated. Expanding ipi about any point on this circle, 
we find tpi cx 5r + iS8, indicating that the nodal circle of "01 represents a vortex 
line. Hence the "01 component forms a vortex ring, within whose core ip2 resides and 
flows azimuthally with the winding number of one. The skyrmion can be viewed as 
a quantized vortex ring in one component, which is filled with the other component 
that also carries the quantized circulation along the ring. This is shown in Fig. 13. 
This configuration closely resembles the "cosmic vortons" which might have been 
formed in the early universe from superconducting cosmic strings 160 . 




Fig. 13. A 3-D skyrmion (left), as a vortex ring supporting a superflow, and a 3-D skyrmion 
(right) with a constant surface density, which was obtained by numerical simulations. The vortex 
core and the ring are shown as the surface of constant atomic density. [Ruostekoski and Anglin, 
Phys. Rev. Lett. 86 3934 (2001), reproduced with permission. Copyright (2001) by the American 
Physical Society.] 

Ruostekoski and Anglin 94 proposed a carefully designed sequence of Rabi tran- 
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sitions that create the skyrmion configuration of Eq. (63). They argued that the 
subsequent dissipative time evolution exhibited unstable dynamics that caused the 
vortex ring component to shrink to zero radius. However, since the general multi- 
component system is characterized by several kinds of two-body scattering lengths, 
whether a stable skyrmion exists or not is a highly nontrivial problem. For homoge- 
neous two-component condensates, Battye et al. found that the energetically stable 
skyrmion can exist under the condition that phase separation occurs, that is, if the 
condition (29) holds 95 . In this case, the two components can strongly repel each 
other and the toroidal filling due to the line component can prevent the vortex 
ring from shrinking to nothing. Consequently, a filled vortex ring, as opposed to an 
empty vortex ring, can be more stable against collapse. Moreover, in the skyrmion 
the filling has one unit of angular momentum about the z-axis, resulting in a 1/r 2 
centrifugal barrier that further hinders the ring from shrinking. 

Savage and Ruostekoski extended the analysis to a harmonically trapped sys- 
tem and found that there are additional instabilities that will play a crucial role in 
experiments 96 . For sufficiently large total number of atoms in the vortex-line com- 
ponent, the nonlinear repulsion between the two components is strong enough to 
inhibit the collapse of the vortex ring, stabilizing the skyrmion against shrinking for 
cylindrically symmetric initial states. However, due to the inhomogeneous potential, 
the skyrmion is still unstable with respect to the drift of the vortex line towards the 
edge of the BEC. Once the vortex line drifts to a low-density region, the nonlinear 
repulsion is no longer strong enough to prevent the shrinking of the skyrmion and 
the collapse occurs as described above. The drift reduces the total angular momen- 
tum of the atoms, indicating an energetic instability. Physically, this results from 
dissipation (e.g., from thermal atoms). Thus, additional physical mechanism, such 
as rotation or optical pinning potentials, will be required for stabilization. 

Rotation of only the vortex-line component gives the simplest example for sta- 
bilization, but this is difficult. Moreover, rotation of both components introduces 
a new instability mechanism. For sufficiently rapid rotations the line component 
again reaches the boundary, resulting in Q < 1. At the same time the ring vor- 
tex accommodated the rotation by twisting and finally breaking at the condensate 
boundary. Nevertheless, for a small range of rotation frequency near 0.085u;_i_, the 
skyrmion was stable against both shrinkage of a vortex ring and drift of a vortex 
line. Then, the stable skyrmion configuration requires a shift of the rotation axis 
to off-center . The breaking of cylindrical symmetry implies a family of degenerate 
skyrmions parametrized by the azimuthal angle. Another method for stabilizing the 
skyrmion is to inhibit the drift towards low-density regions by creating a positive 
density gradient around the vortex line. This can be implemented with a blue- 
dctuned Gaussian laser beam along the z axis, providing a cylindrically symmetric 
repulsive Gaussian dipole potential perpendicular to the z axis. Under these condi- 
tions, we expect that atomic BECs, which have been imprinted with a topological 
configuration, will relax via dissipation to a stable skyrmion configuration. 
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4.5. Dynamic properties of vortices 

The dynamic properties of vortices in two-component BECs are also expected to be 
extremely rich because of the various kinds of interactions between the two super- 
fluids. A remarkable feature of this system arises from "buoyancy" 161,162 ' 85 , which 
refers to a net mean-held force on a constituent of a multicomponcnt condensate, 
the force being dependent on the various intraspecies and interspecies interaction 
parameters. 

As shown in Sec. 4.1, the unstable dynamics of the vortex was observed, depend- 
ing on which component has a vortex in the 87 Rb condensate 11 . This phenomenon is 
understood from linear stability analysis of the vortex state based on the mean-held 
approach 82,83,84 . We use the following notation for the states: (1,0) for the state 
with the vortex in |1) and (0, 1) for the state with the vortex in |2). In the JILA 
experiment, two components had the same number of particles (Ni = N 2 = N) but, 
in general, one could consider any ratio between the populations of the different hy- 
perhne states. For 87 Rb atoms, the relation of the coupling constants g\ > gn > g 2 
(ai : a.12 : a 2 = 1-03 : 1 : 0.97) favors a configuration that has its hrst component 
spread over the largest part of the space, as shown in the previous discussion of 
equilibrium. Numerical simulations show that for equal populations N± = N 2 = N , 
and with the above coupling constants, the stationary states (1,0) is stable while 
the other state (0, 1) is unstable 82,83,84 . 

The origin of the instability of the state (0, 1) is purely dynamical, which can be 
understood within the framework of the mean-held theory for the two-component 
system without dissipation. Actually, the instability does not lead to expulsion of 
the vortex from the condensate, but instead to periodic transfer of the vortex from 
one component to the other. To study the vortex stability, one starts from the time- 
dependent, coupled Gross-Pitaevskii equations for the condensate wave functions of 
each species: 

h 2 V 2 



thr w 



ih 



dt 



2m 
h 2 V 2 
2m 



+ ^i+ 5 i|*i| 2 + ff i 2 |*2| 2 
+ V 2 + .g 2 |* 2 | 2 + Si 2 |*i| 2 



*i, (64) 
* 2 - (65) 



These equations are the particular case of Eqs. (31) and (32) in which the drive is 
turned off (Hi, wr, A = 0). Equations (64) and (65) conserve the number of particles 
in each hyperhne state. However, the angular momentum of each component is no 
longer a conserved quantity. Instead, the conserved quantity is the total angular 
momentum (L z ) = J d 3 r(^lL z ^/i + ^ 2 L Z ^ 2 ). As in the JILA experiments, both 
potentials are assumed to be spherically symmetric and have the form Vi(r) = 
^(r) = \mLo\ir 2 + z 2 ) . For stationary configurations in which each component has 
a wcll-dchncd value of the angular momentum, the time and angular dependences 
of the wave functions are factored out as 



*i(r, z, 0) = e-^ h e^ e Mr, z) (i = 1, 2). 



(66) 
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Here, we focus on the two axisymmetric configurations with vorticity (71,72) = 
(1,0) and (0,1), which correspond to the single vortex states for the |1) and |2) 
components shown in Fig. 3. 

Linear stability analysis of the two axisymmetric vortex states leads to the fol- 
lowing results. Among the normal modes of the (1,0) case, there is a negative 
eigenvalue, which means that there is a path in the configuration space along which 
the energy decreases. This path, which is associated with the mode of the dipolc 
type, belongs to a perturbation that causes the initial displacement of the vortex 
from the trap center. However, this drift instability of the vortex is affected by dis- 
sipation; in particular, without dissipation, the configuration is dynamically stable. 
But, it can be energetically stabilized by turning on rotation as discussed in Sec. 
4.2. In the (0, 1) case, on the other hand, there are normal modes with complex fre- 
quencies. The unstable modes have a similar shape of the negative-eigenvalue modes 
of the (1, 0) case, which means that the perturbations also displace the vortex from 
the center. The imaginary part of the eigenvalues implies that vortices with unit 
charge in |2) are dynamically unstable under a generic perturbation of the initial 
configuration and this instability grows without dissipation. This result is consistent 
with the JILA experiments in Fig. 2, in which a vortex in the |2) component was 
unstable and collapsed. 

Numerical simulations of the vortex dynamics for large perturbations were done 
by Garcia-Ripoll and Perez-Garcia 82 . The results showed that the linearly stable 
state (1,0) is robust and survives under a wide range of perturbations; the vortex 
only shows a precession around the center. In contrast, the unstable configuration 
(0, 1) gives rise to recurrent dynamics: the first component and the vortex oscil- 
late synchronously, where the vortex core in |2) pins the density peak of |1). These 
oscillations grow in amplitude until the vortex spirals out. The displaced vortex 
carries less than the total angular momentum, and so the vortex precesses around 
the center. Because of this and the conservation of the total angular momentum, 
the vortex is transferred from |2) to |1). Although the dynamics is not completely 
periodic, this mechanism exhibits some recurrence and the vortex eventually returns 
to 1 2). The preceding behavior persists even for strong perturbations in a 2-D con- 
densate. However, for large perturbations of a 3-D condensate, the dynamics may 
be chaotic. 

The preceding results are valid when iVi = N 2 - For an arbitrary ratio of the 
populations N1/N2 and arbitrary values of the coupling constants gi, gi and 
P12, the dynamical stability conditions can be estimated analytically 83 ' 84 . Un- 
der the two- mode approximation 83 , where the wave function is approximated as 

~ ^(t)^ g (r) + bi(t)^! e (r) with the cigenfunctions of the d-dimensional harmonic 
oscillator * s (r) = (l/7r) d / 2 e - r2 / 2 and * e (r) = (2/ ' d-nf^re'^ '/ 2 e lt> , the configura- 
tion (1, 0) is stable if 




(67) 
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For 87 Rb, the inequality (67) is always satisfied because the right-hand side is neg- 
ative, which proves that the configuration with a vortex in |1) is always linearly 
stable. Note that the stability properties do not depend on the total number of 
particles but only on the ratio between the populations. The stability condition of 
the configuration (0, 1) can be obtained in a similar manner as 



This inequality fails for a certain range of N1/N2. For the case of 87 Rb the unstable 
range is 0.68 < N 2 /Ni < 1.38, which means that certain choices of the population 
imbalance allow stabilization of the vortex in |2). More detailed calculations were 
done by Skyrbin 84 who derived the stability condition as 



which gives an unstable range of 0.63 < N 2 /Ni < 1.69, which is slightly different 
from Eq. (68). These results predict the possibility to dynamically stabilize the 
axisymmctric vortex states for various multiple-condensate systems. 

The effect of the intercomponent interaction was also demonstrated in a pre- 
cession of an off-centered vortex core 161 . This is a consequence of a Magnus effect; 
that is, an applied force on a rotating cylinder in a fluid causes pressure imbalances 
at the cylinder surface that makes the cylinder drift orthogonal to the force. The 
experiment observed that the precession frequency of the filled core is slower than 
that of the empty core 161 . The slower precession of filled cores can be understood 
in terms of the buoyancy effect. Because of its slightly smaller scattering length, 
the 1 2) component has negative buoyancy with respect to the |1) component, and 
consequently tends to sink inward towards the center of the condensate. With in- 
creasing amounts of the |2) component in the core, the inward force on the core 
begins to counteract the outward buoyancy of the vortex velocity field, resulting in 
a reduced precession velocity. It is predicted that with a filling component of suffi- 
ciently negative buoyancy in the core, the core precession may stop or even precess 
in the direction opposite to the direction of the fluid flow 162 . 

4.6. Quantized vortices in spin-1 Bose-Einstein condensates 

Another multicomponent Bose-condensed system is a "spinor BEC" characterized 
by the hyperfine spin F = 1. An optical trap removes the restriction of the confin- 
able hyperfine spin states, realizing new multicomponent supcrfluids with internal 
degrees of freedom 71 ~~ 76 . In this section, we review briefly the study of quantized 
vortices in such a BEC with |F| = 1 hyperfine spin states |F=1,tof = ±1,0). 

To describe the F = 1 spinor BEC, a generalized GP mean-field model was 
introduced by Ohmi and Machida 97 and Ho 98 . With the mean-field approximation, 




(68) 




(69) 
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the total Hamiltonian of the system reduces to an energy functional similar to Eq. 
(14): 

**] = / [e ** (-^r + V ^ - **+|E 

(70) 



f E E *i ^(^iA)**^ 



Here we consider the system in a rotating frame with f2 = fiz and the subscripts 
represent a = (x,y,z) and i,j,k,l — (0, ±1). The condenate wave functions are 
characterized by the hyperfme sublevels rap = 1,0,-1 as * = (\&i, vf^, *_i). If 
one of the component ^ is fixed at zero, Eq. (70) reduces to the energy functional 
of the two-component case Eq. (14); the crucial difference from the two-component 
problem is the presence of the g s -term. The coupling constants g n and g s are written 
as 

4irh 2 2a 2 + a a Airh 2 a 2 - a , 

9n = 5 , 9s = z ' 71 

m 3 to 3 

where an and a 2 are the scattering lengths for two colliding atoms with total spin 
angular momentum F tot = and 2. The spin operators F a (a — x,y,z) can be 
expressed in matrix form as 





F z = . (72) 



Note that the term proportional to g s in Eq. (70) describes the spin exchange 
interaction, which is ferromagnetic for g s < and antiferromagnetic for g s > 0. The 
spinor BEC of F — 1 23 Na atoms studied by the MIT group was demonstrated to 
have an antiferromagnetic interaction 71 , while the F = 1 87 Rb condensate has a 
ferromagnetic one 75 ' 163 . The sign of g s plays a crucial role in determining the stable 
vortex structure. The time-independent coupled GP equations are derived by the 
requirement that the energy functional (70) be extremal: 

^ ' i 

+.9s*-i(*o) 2 = (73) 

(-^- + ^( r ) - *o + (g n + s.)(E *?*i)*o - ,g s *S*o*o 

+2. 9s *** 1 *_ 1 - Mo * , (74) 
V(r) - QL^ *_! + ( 5 „ + ff8 )(^ - 2 5s *** 1 vl/_ 1 

' i 

+ffa**(*o) 2 = (75) 



ft 2 V 2 
2to 
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where the three chemical potentials in (i = 0, ±1) play the role of Lagrange 
multipliers. They are determined so as to conserve the total particle number 
N = J dr(|*i(r)| 2 + |* (r)| 2 + l*-i(r)| 2 ) and magnetization M = Jdr(|*i(r)| 2 - 
|^_i(r)| 2 ). The angular momentum (L z ) = J2i I dv^iL^i should also be consid- 
ered as a conserved quantity. 

4.6.1. Coreless vortex states 

First consider the structure of axisymmetric vortex states. The detailed analysis 
was made by Isoshima and Machida 103 . In this structure, the vortex configuration 
is classified by the combination of the winding number ji of the condensate wave 
function tyj — \f~f>j(rj exp[i(aj +"/j0)] with the cylindrical coordinate (r, 9, z), where 
the homogenity of the wave functions along the z-axis is assumed, and oij is the 
overall phase of the j-th component. The phase factors of are determined such 
that the energy of the g s term in the integrand of Eq. (70) is minimized; this term 
reads 

a ijkl 

+ (pi-P-i) 2 }- (76) 
To minimize E s , ctj and should satisfy 

2a = a>i + q_i + nir, (77) 
270 = 71+7-1, (78) 

where n is an integer, and the odd (even) n corresponds to the antifcrromagnetic 
(ferromagnetic) interaction. The global phase ctj has no effect on the vortex struc- 
ture, and will therefore be set to zero. The possible combination of 7, satisfying 
Eq. (78) gives this system a characteristic vortex structure. If the values of jj are 
restricted to being -fj < 1, we have (71, 70, 7-i) = (1,0,-1), (1, 1, 1) and (1, 1/2, 0). 
In the case of (1, 1/2, 0), the value 1/2 is not allowed in the axisymmetric system, so 
that the ^ -component vanishes; this vortex state is called the "Alice state" (half- 
quantum vortex state) 79 ' 103 . Therefore, there is one-to-one correspondence between 
the Alice state and the axisymmetric vortex state in two-component BECs discussed 
in Sec. 4.2 (i.e., *i(*_i) — > i^ife) in Fig. 4). The structure and the stability of 
these vortex states for various ft and the total magnetization M were thoroughly 
investigated in Ref. 103. 

Also, the winding number can exceed unity. Mizushima et al. discussed the 
(0, 1, 2) vortex shown in Fig. 14 and showed that the ferromagnetic interaction sup- 
ports the thermodynamic stability of this structure 104 . The central region of the 
harmonic trap is occupied by the non-rotating ^1 component. The ^0 component 
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with 70 = 1 is pushed outward, while the component with 7_i = 2 occupies 
the outermost region. The resulting total density is nonsingular. This configuration 
is favorable for the ferromagnetic case, because it is favorable to spatial phase sep- 
aration; the presence of vortices with different jj effectively cause phase separation 
in the radial direction as shown in Fig. 14(a) 




(a) (b) 

Fig. 14. A (0, 1, 2) vortex in a ferromagnetic spinor BEC. The density profile of the condensates 
(a) and the density map of the / z -vector (b) for the (0, 1, 2) vortex at CI/lu± = 0.35 and M/N=0.21. 
In (a), the bold curve is the total density px = 2Zj l^j| 2 an d the thin curves show the density of 
each component |*j| 2 (j = 1, 0, -1). [Mizushima, Machida, and Kita, Phys. Rev. Lett. 89 030401-1 
(2002), reproduced with permission. Copyright (2002) by the American Physical Society.] 



Putting this coreless vortex into the pseudospin representation reveals that it has 
the structure of a spin texture. The wave functions of Fig. 14(a) can be parametrized 
as 

/ 2 0(r) 




vW(r) 



V 



e ie^ sin m cos en | , (7!)) 



where the bending angle (3(r) has the range < (3(r) < n. The spin direction is 
denoted by the /-vector and is given as l(r) = z cos (3(r) + sin/3(r)(cos Ox + sin Oy) 
[the contour density of l z is shown in Fig. 14(b)], where [3 varies from /3(0) = 
to (3(R) = ir/2 (= 7r) for a Mermin-Ho (Anderson- Toulouse) vortex (R is the 
outer boundary of the cloud). In other words, l z = cos (3(r) varies from l z (0) = 1 
to l z (R) = (—1), where i-vector points up at the center and points outward 
(downward) at the circumference for the Mermin-Ho (Anderson- Toulouse) vortex. 
In the superfluid 3 He-A phase, the stability of these vortices is guaranteed by the 
constraint that the /-vector is perpendicular to the vessel wall 2 . There is no such 
boundary condition in the trapped BEC system. In this case, the (0,1,2) vortex 
of Fig. 14 can be interpreted as a vortex with an intermediate boundary condition 
(tt/2 < (3(R) < 7r); wc can control the value of (5{R) from the Mermin-Ho condition 
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to the Anderson- Toulouse one by merely changing the total magnetization M 104 . 
In another numerical study, Martikainen et al. analyzed the coreless vortex state 
as a function of rotation frequency without fixing the total magnetization, and 
found that (3{R) increases with f2 and the upper value of (3(R) is 37r/4, above 
which additional vortices nucleate. This implies that a Anderson- Toulouse vortex 
can never be the ground state of the system. 

Full 2-D numerical simulations of the coupled GP equations (73)-(75) also 
reveal non-axisymmctric structures of vortices, whose stability depends sensi- 
tively on the spin-exchange interaction g s , the external rotation O and the total 
magnetization 99,104 ' 105 . For example, the axisymmetric (1,1,1) vortex is always 
unstable, leading to the displacement of the vortex cores from the center. In the 
antiferromagnetic regime, the vortex cores are displaced such that the condensates 
has two singularities, through an overlap of the and ^-i components or it has 
three singularities to form a triangle configuration. The breaking of the axisym- 
metry mainly leads to a decrease in the overlapping area of the condensate wave 
functions. In comparison with axisymmetric states, these non-axisymmetric vortices 
have the advantage that they can easily adjust to changes in Q. As fi increases, the 
two or three separate singularities adjust their interdistances from the center and 
change the value of L z to gain energy due to the term — QL Z . In this sense, the 
non-axisymmetric vortices are more flexible for a change in compared with the 
axisymmetric ones. 



4.6.2. Monopoles 

Another interesting topological excitation that occurs in the spin-1 BEC system is 
a singular pointlike defect, called a "monopole" due to its similarity to magnetic 
monopoles in particle physics 164 . A monopole is characterized by a unit vector that 
is radial with respect to a unique central point (i.e., the "hedgehog" structure). 
Stoof et al. 113 proposed monopole excitations in an antiferromagnetic spin-1 BEC. 
The monopole is characterized by the spinor 

/ \ H^M ( -m x {r)+im y {v)\ 

*o = f-^fe^) V2m z (r) , (80) 

V m x (r)+im v (r) J 

where 9 is the overall phase. The form of Eq. (80) originates from the fact the 
antiferromagnetic interaction implies zero average spin |(F)| = 'ipF^/px = 0. The 
spherically symmetric monopole is described by the radial unit vector m(r) = r/r. 
Minimization of the energy of the symmetric configuration yields 9 = and the den- 
sity singularity at the origin /Ot(0) = 0. Then, the spinor component Wo resembles 
a dark soliton and form perfectly overlapping, straight singly-quantized vortex 
lines with opposite circulation, the latter being perpendicular to the phase kink 
plane. Stoof et al. demonstrated that this particular spin texture is a consequence 
of the unit winding number and the minimization of the gradient energy 113 . This 
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configuration could be created through the phase imprinting technique achieved by 



Fig. 15. A stable half-quantum vortex ring (Alice ring). To form it, the energy of a spherically 
symmetric monopole was minimized by continuously deforming the field configuration. The cal- 
culation was done for a condensate trapped in a spherically symmetric trap with the frequency 
uj = 2tt X 10 Hz and the total particle number TV = 4 X 10 6 in a weakly antiferromagnetic regime 
9s 1 9n = 0.04. The left figure shows the asymptotic distribution of the spin quantization axis m(r); 
for r 3> £a it forms the radial hedgehog. For visualization purposes, the unoriented m(r) field is 
drawn by cones. The right figure shows the constant surface density for |<E'i(r)| 2 (light grey) and 

for | i (r) | 2 (dark grey), where the monopole core is deformed for r < £ a with the two line 

vortices separating. [Ruostekoski and Anglin, Phys. Rev. Lett. 91 190402-1 (2003), reproduced 
with permission. Copyright (2003) by the American Physical Society] 

The theoretical study by Ruostekoski and Anglin 116 demonstrated that the 
spherically symmetric monopole is stable only in the strongly antiferromagnetic 
regime. This is due to the following two length scales 



They represent the length scales over which pt and | (F) | , respectively, tend to their 
bulk value when subjected to localized perturbations. In the weakly antiferromag- 
netic regime, the wavelength (~ £ a ) at which the antiferromagnetic constraint may 
be violated is much larger than that (~ £ s ) at which the total density constraint 
fails. Then, the monopole core extends to a larger size, which is favorable for having 
nonzero average spin instead of a density zero. The monopole point defect may be 
continuously deformed into a circular line defect without changing the topological 
invariant. This is done by punching a hole in the spherically symmetric core. To 
keep $ single-valued on the disc bounded by the ring, the macroscopic phase 9 
must change by n around any loop that links the defect circle, while there is also 
a 7r-disclination in m on the disc. This structure is identified as a half-quantum 
vortex line, forming a closed circular ring, called an "Alice ring" [see Fig. 15]. To 
keep \& single-valued on the ring itself, one can either have px vanish there, or have 
|(F) | = 1 on the ring instead. By decreasing the ratio g s /g n , the spherically sym- 
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metric monopole core with pr(0) = becomes unstable against formation of an 
energetically stable Alice ring with |(F)| = 1 core. According to a linear stability 
analysis of this case, the radially symmetric monopole is unstable in a homogeneous 
condensate for g s /g n < 0.17, while according to the numerical simulation for a 
trapped condensate the monopole is unstable for g s /g n < 0.16, using an experimen- 
tally feasible set of parameters for 23 Na (uj± — oj z — 2tt x 10 Hz and N — Ax 10 6 ). 

The creation of the monopole is not limited to the antiferromagnetic spinor con- 
densates. Martikaincn et al. studied in detail the two-component monopole charac- 
terized as 

/ \ IVJr) / -m x {r) + im y (r)\ 
: „ Wf^ V2m z (r) J. (82) 

To ensure the stability of this texture against phase separation, the spin-1 conden- 
sate must have ferromagnetic interactions, which makes the 87 Rb spinor conden- 
sate a potential candidate. In other words, contrary to a previous expectation, the 
preparation of a monopole is not restricted to the antiferromagnetic 23 Na conden- 
sate because one can use spinors that do not have the order-parameter space of the 
ground state. 

5. Vortex states in rapidly rotating two-component BECs 

Rapid rotation of a BEC generates a large number of quantized vortices, with the 
result that the condensate can mimic a rigid body rotation despite the irrotational 
nature of superfluidity 1 . It is well known that the most stable configuration of an 
assembly of quantized vortices in a rotating single-component superfluid ( 4 He-II 
or atomic gas BECs) is a triangular lattice; also it is known that this is analo- 
gous to an Abrikosov lattice of magnetic fluxes in a type-II superconductor. At 
present, most experiments on vortex lattices in atomic-gas BECs have been done in 
single-component systems. Thus, a natural question arises: what happens in rapidly 
rotating two-component BECs? The vortex lattices in such systems are more com- 
plicated than those in single component condensates because of the intercomponent 
interactions. This section reviews the studies of vortex states in two-component 
BECs with a large number of vortices. We also describe a recent experiment 165 that 
found interlaced square vortex lattices in two-component BECs. 

First, we show results obtained through analysis based on the "mean-field quan- 
tum Hall regime" 87 . This regime will be valid when the BEC in a harmonic trapping 
potential is rotated with angular frequency ft close to the radial trapping frequency 
ClI± 43,166_ This method has the great advantage that one can analytically describe 
the vortex structure, although the method is effective only near the limit f2 ~ u)j_. 
Next, we show the results of the direct numerical calculation of the coupled GP 
equations (21) and (22) 88 . They reveal a rich variety of vortex structures that have 
eluded analytic treatment. Depending on the strength of the intercomponent inter- 
action and the rotation frequency, vortex lattices can have the form of triangular 
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lattices, square lattices, double-core vortex lattices, and "serpentine" vortex sheets. 
These properties can be understood by using the concept of "pseudospin" , in which 
the vortex lattices are interpreted as a lattice of skyrmion. Finally, we discuss the 
effect of internal coherent coupling on the lattice structure, caused by the anisotropy 
of the resulting vortex molecules (meron pairs). 

Similar discussions can also be done for spin-1 BECs. Kita et al. studied vortex- 
lattice structures of antiferromagnetic spin-1 BECs using the phenomenological 
Ginzburg-Landau equation 106 , which is formally similar to the mean-field quan- 
tum Hall approach. Their study revealed that the conventional Abrikosov lattice 
with hard vortex cores is unstable and the vortex cores shift their locations. The 
system has many metastable configuration of vortices, depending sensitively on the 
ratio (72/30 in Eq. (70). The vortices are characterized by distribution of magne- 
tization and the difference in the number of circulation of quanta per unit cell, 
all of which makes the characterization of the vortex lattice structure complicated. 
Further discussions on the phase diagram of the ground state in rotating spin-1 
bosons are reported in Reijnders et al. 109 . This paper also describes a study of 
vortex-lattice structures with the mean-field quantum Hall approach and the exact 
diagonalization study of the quantum liquid phase. In addition, Mizushima et al. 
did numerical simulations of the GP equations (73)- (75) for a ferromagnetic spinor 
BEC in a fast-rotating regime and also in an external magnetic field 111 . As a re- 
sult, they could constuct a phase diagram of stable configurations on a plane of the 
rotation frequency and the total magnetization. 

5.1. Mean-field quantum Hall approach 

The "mean-field quantum Hall regime" is the regime obtained in a rapidly rotating 
limit, where mean-field theory remains valid. Thus, each component can still be 
characterized by a condensate wave function (i = 1, 2). The angular momentum 
of the system is so high that \&i is made up of the orbitals in the lowest Landau 
level (LLL) in the plane perpendicular to the rotation axis. Ho showed that this 
regime will emerge in a 3-D Bose gas at sufficiently high angular momenta 166 , and 
recently the JILA group has achieved this regime experimentally 43 . (See also the 
recent detailed arguments in Ref.167-170 for the use of LLL approximation and 
the structure of a vortex lattice.) In this regime, the wave function acquires an 
analytic structure that allows exact evaluation of the energy of a vortex lattice. 
As a result, it is possible to evaluate a wide range of lattice structures in rapidly 
rotating condensates, a feat that would be impractical for numerical calculations 
because of the time and the accuracy required. 

An extension of this argument to the two-component BECs was made by Mueller 
and Ho for the following situation. If each component has an equal number of 
bosons, and if the trapping potentials of the two components are identical, the two 
components will have the same size and the same density of vortices. In this case, 
one expects that each component will contain identical vortex lattices, with one 
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lattice displaced relative to the other because of the presence of the intercomponent 
repulsive interaction. We shall refer to such a configuration of vortex lattices as 
staggered vortex lattices. While in the current experiment, the coupling constants 
satisfy <?i ~ 172 ~ 512, considerable insights are gained by studying vortex lattices 
for a wider range of coupling constants. Mueller and Ho considered the case g\ ~ 
.92 7^ .912 and investigated a wide range of vortex lattice structures as a function of 
8 = 912/ \j 9i92, and we summarize their result here. 

The condensate wave functions ^1 and ^2 of a rotating two-component BEC 
are determined by minimizing the free energy F = E — /i X N\ — where E is the 

total energy functional (14), fii (i — 1, 2) are the chemical potentials that determine 
the number of bosons iV, in each component. For simplicity, we assume the same 
cigar-shaped trapping potential V — m{bj\r 2 + uj 2 z 2 )/2 with ui± >• u) z for each 
component. As discuss in Ref.166, a slow variation of the trapping potential along 
the z-axis allows one to apply a Thomas-Fermi approximation for the z dependence 
of 'J, and write F as 



J ^1=1,2 



1 (K. 



, V — mwj_z x r + (wj_ — 0,)L Z — Hi(z) 
2m \ i 



+ ^l|*l| 4 + ^2|*2| 4 +5l2|*l| 2 |* 2 | 2 > (M) 



with i-ii(z) = in — mLo 2 z z 2 /2. As z is treated as a parameter, it is convenient to write 
*i = \/Pi{z)ipi(r, z) with J \ipi(r, z)\ 2 d 2 r = 1; hence, the normalization condition 
J d 3 r\^>i\ 2 = Ni is now written as J dzpi(z) = A^. 

The first term of the integrand in Eq. (83) /il = (— ifiV — muj±z x r) 2 /2m is 
identical to the Hamiltonian that describes the motion of an electron in a magnetic 
field Bz in the symmetric gauge (i.e., (— iSV — eA/c) 2 /2m with a vector potential 
A = Bz x r/2). The eigenstates of are referred to as the Landau levels and 
described in polar coordinates (r, 9) as 



u 



n ,m 



/ \ \m—n\ / 2 \ 

(r) - O P *(«-™)e-r 2 /4&L { J_ \ j-\rn-n\ / J_ \ 

(r)~U nm e \b ho J (rn+n-\m-n\)/2y 2 b 2 )' [ ' 



where 6ho = y fr/mui±, C nm is the normalization constant, and L 1 / n {x) the associ- 
ated Laguerre polynomial. Since u n>m are also eigenstates of L z with the eigenvalue 
h(m — n), the eigenenergy of hi, — (ui± — Cl)L z reads 

c nm = h(u>± + fl)n + h(u>j_ — fl)m + huj±, (85) 

where n is the index of the Landau level, and m is the index that classifies the 
degenerate states in the n-th Landau level. This degeneracy with respect to m is 
lifted by the interaction. When approaches ui± 7 the condensate density becomes 
so dilute because of the centrifugal effect that the mean-field interaction energy 
decreases. If the interaction energy becomes less than the energy spacing of the 
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Landau level, then the wave functions ipi can be written in terms of the orbitals 
u o,m( r ) hi the LLL in the xy-pl&ne, 



«0,m0) 



1 



\[2jrm\ 

The energy functional then becomes 



x + iy 
bho 



dz 



E 

i=l,2 



h(u± — ft) 



Ko 



fj,i(z) + huj_ 



(86) 



Pi{z) 



+ j d2r \ E + 9i2Pi P 2m 2 m 2 } (s?) 

where (r 2 ); = fr 2 \ip t \ 2 d 2 r. 

We next introduce a variational wave function that represents a vortex lattice. 
Generally, wave functions in the LLL can be written as 0(r) = f(z)e~ r / 2b h°, where 
f(z) is an arbitrary polynomial of z = x + iy (but not of z*). The form of a vortex 
lattice can be described by 

2 /2fcL 



<f>(r) = C]J(z-b a )e- 



(88) 



where C is a normalization constant and {b a } are the zeros of <j>, which are the 
positions of vortices. In Eq. (88), we assume perfect periodicity of the lattice. If 
the zeros form an infinite lattice with unit cell size v c , the density \<fi\ 2 can then be 
written as a product of a Gaussian and a function which is periodic under lattice 
translation: 



I 2 = e-'' 2 /- : 



g(r), g{r)=g{r + TL), 



(89) 



where R = riiBx + n 2 B 2 , rij are integers, and B l7 B 2 are the basis vectors of 
the lattice. The width a includes the contribution of the number of vortices of the 
system, which is determined below. 

It is convenient to introduce the complex basis vectors (61,62) defined by 6, = 
(x + iy ) ■ Bi ; then the unit cell size is v c = i (61 6 2 — 6* 62 ) /2 . If the lattice is oriented 
so that 61 is real, we can parametrize Bi = 61X and B 2 = b\(uic + vy), and obtain 
62 = bi(u + iv) and v c — b 2 v. Except for the normalization constant, the form of 
Eq. (88) is similar to the Weierstrass's cr-function 171 



a(z) = z [] 



1 - 



a 



cxp 



S 1 > ■ 



+ 



2m 



(90) 



where the term with m = n = is omitted from the product and Q mn — 2muj 1 + 
2nui2 with the half-periods w\ and u> 2 - By using a relation between the Weierstrass's 
cr-function and the Jacobi ^-function 171 



0(C, 



E i-^e 



i7rT(n+l/2) 2 +27riC(n+l/2) 



(91) 



l 



l 



2, 2008 15:1 WSPC/INSTRUCTION FILE 2vorrev 



48 K. Kasamatsu, M Tsubota and M. Veda 

we can write f{z) = 6((,r)h(Q, where £ = (x + iy)/b\ = x + iy, r = u + iv = 6 2 /&i, 
and h(Q is an entire function without zeros. To ensure the normalization of <f>, the 
function h(Q must be of the form h{Q = exp(ciC + c 2 £ 2 ). 

The density of the system is then |0(r)| 2 = \9((,T)\ 2 \e c ^ +c ^ 2 \ 2 e- r2 ^ . For a 
vortex lattice with inversion symmetry about the origin r = 0, we have c\ = 0. 
In addition, if the shape of the condensate is cylindrically symmetric, we have 
c 2 = n/(2v c ), which gives |0(r)| 2 = \9(£, r)| 2 e - r2 / <j2 with 

Following Eq. (89), we can write g(r) as 

g(T) = \0(t,T)e-"v'' v °\ 2 . (93) 
It follows from the quasi-periodic properties of the Jacobi thcta-function 

0(C + l,r)=fl(C,r), (94) 

#(C +t,t) = -e-^+^OiC, t) (95) 

that g(r) = g(r + H). The periodicity of <7(r) implies g(r) — v~ x ^2 K ffKe lK r > where 
{K} are the reciprocal lattice vectors. K = miKi + m 2 K 2 , and Kj are the basis 
vector of the reciprocal lattice, Ki = (27r/u c )B 2 x z, K 2 = (2n/v c )z x Bi, and 

Vc K 2 = ((vmi) 2 + (m 2 - umi) 2 ) . (96) 

v 

The explicit form of gn is calculated in a straightforward way. We obtain 

r ^|2 _ V ' e i7ru[m 2 (m2+l)— mi(mi+l)] e -7TO[m2(m 2 +l)-mi(mi+l)+l/2] 

mi ,7712 

x e 2i7rx(m 2 — mi ) e —27ry (mi +7712 + 1) 

= E( _1 ) me2ff<m * e_,r '' m2/2Lm > ( 9? ) 

m 

where 

£ ^_ ^ ^ ^ g?7r(m+m')^ ^(i-Trum — 2iTy — -nvm' /2}m' (98) 

m' 

The factor (1 - e ^(™+™'))/2 is added in Eq. (98) because m + m' must be odd. By 
applying the Poisson summation formula 

+ 00 +00 p + OC 

E /(")= E / dxf(x)e^ k * (99) 

to Eq. (98), we have 

7- _ ' [ -Tt(2k+um+2iy) 2 /2v , /_ i \ra+l -7r(2fc+l+«m+2iy) 2 /2« 

m ~V2^v^ ( j J 

_ 1 y^ |-_ 1 ^(m+l)fc e -7 r (fc+Mm+2i7j) 2 /2u^ ^qq-j 
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We thus have 

|e(»)l 2 



m,n 
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m+n+mn 
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K 



2ny /v c 



where r = xx + yy, K = (27rmx — 2n(n + um)/vy)/b\, and 

gK = (_l^m+n+mn e -v c \K\ 2 /8ir 



(101) 



(102) 



If gi = g 2 , the two components with equal particle numbers and trapping poten- 
tials have identical vortex lattices. A sufficiently small difference in gi — y 2 should 
not change this structure. If the intercomponent interaction is repulsive </i 2 > 0, 
the mean-field interaction energy is minimized by interlacing the two lattices; if the 
vortex lattice of one component is deformed, the other has to follow to keep the 
interaction energy minimal. Under these assumptions, the normalized condensates 
ipi and ip2 in Eq. (87) would be represented by densities of the form 



i^^E^-e- 

7T(T Z z ' 

9K 



K 



E K <SK' e - 2K ' 2 / 4 ' 



(103) 
(104) 
(105) 



The wave function is described by variational parameters Pi(z), a 2 , the basis vectors 
Bi (which determine the unit cell size v c ), and the relative displacement ro between 
the two lattice. 

By integrating Eqs. (103) and (104), one would see that terms up to those of 
relative order v c /a 2 , the radius of the condensate is (r 2 ) x = (r 2 ) 2 = a 2 . Defining 
the quantities / and I12 as / <i 2 r|<I>i| 4 = I/(na 2 ) and J d 2 r|$i| 2 |$ 2 | 2 = Ii2/(^cr 2 ), 
we have 



/4E2KSK< e - CT2 ' K+K '' 2/8 , 



K,K' 



/i2=E5K.9K'e- 4K '- r " e - CT2|K+K ' |2/8 , 



(106) 
(107) 



K 



and the energy functional F takes the form 



F 



/4 



fi(z) — hu)± — h(ui± — fi) 
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(108) 



2, 2008 15:1 WSPC/INSTRUCTION FILE 2vorrev 



K 



9k 

go 



50 K. Kasamatsu, M Tsubota and M. Veda 

Since we assume here a large number of vortices, the size of the condensate is 
much larger than the unit cell and thus ira 2 /v c » 1. We can therefore ignore all 
K + K' ^ terms in Eqs. (106) and (107), since <r 2 K 2 > ir<r 2 /v c . We then have 

2 

cosK • r , (109) 

K y0 

where <7k is given by Eq. (102) and the sum of K is over the integers mi and 
m 2 . Since the expressions of I and ii 2 in Eq. (109) are independent of er 2 , the 
minimization of F in Eq. (108) with respect to a 2 and pi is greatly simplified. The 
optimized a 2 and p — (p\,p2) are given by 

2 2 n(z) - hw ± 

= 6ho 3^-f7)' (U0) 
2a 2 

PW = ^3-M-^±]G- 1 -1 ) (111) 

where 

G =U'/,/!?)' -(!)• (112 » 

and the corresponding free energy F is given by 

F = - J dz^ \fi(z) - fiwi] 1 • p(s), (113) 

It is clear from Eqs. (110) through (113) that the solution for the case where g\ — 
92 -C | ,9i + 52 1 is very close to that of g\ = g^. The lattice shape (parameterized by 
r , u, and v) enters the grand potential only through the factor 1 • G _1 • 1. When 
.9i = 92 this term is inversely proportional to J = I + 5I\2 (<$ = S12/.91), an d the 
most favorable lattice is the one that minimizes J. 

In the minimization it is convenient to measure lengths in unit of the the basis 
vector Bi = bix and also to write the complex representation of B2 and ro as 
r = u + iv = \T\e lri and vq = a + br, respectively. The phase diagram of the vortex 
lattice as a function of the ratio S = g\ilg is shown in Figs. 16. The main features 
are 

• (a) 6 < 0: In this region the vortices of the two components coincide with 
each other (a = b = 0) to form a triangular lattice (t = e lTr / 3 ). 

• (b) < S < 0.172: The vortex lattice in each component remains triangular. 
However, the relative displacement r between the two lattices undergoes a 
first-order change so that one lattice is displaced from the center of the tri- 
angle of the other component (a = b = 1/3). The lattice type (characterized 
by t = e l7r / 3 ) remains constant within this interval. 

• (c) 0.172 < S < 0.373: At 5 = 0.172, r moves from the center of the triangle 
(i.e. half of the unit cell) to the center of the (rhombic) unit cell (a = b = 
1/2). The angle rj abruptly changes from 60° to 67.95° at S = 0.172, and 
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Fig. 16. Vortex lattice phases in rotating two-component BECs. Black and grey dots represent 
vortices of each of the two BECs. The panels (a) through (e) show the vortex structure in each 
of the phases described in the text. The final panel depicts the lattice types; the black and grey 
dots, respectively, occupy positions in the complex plane {m + nr} and {(a + m) + (6 + n)r}, 
where m and n are integers. All minimal-energy configurations have a = b. [Mueller and Ho, Phys. 
Rev. Lett. 88 180403-3 (2002), reproduced with permission. Copyright (2002) by the American 
Physical Society.] 

increases continuously to 90° as 5 increases to 0.372. As a result, the lattice 
type is no longer fixed and the unit cell is a rhombus. The modulus of t, 
however, remains fixed across this region. 

• (d) 0.373 < 5 < 0.926: The two lattices are "mode-locked" into a centered 
square structure throughout the entire interval (r = i, a = b = 1/2). 

• (e) 0.926 < 5: The lattice type again varies continuously with interaction 
S. Each component's vortex lattice has a rectangular unit cell (rj — tt/2) 
whose aspect ratio |r| increases with S. Both 87 Rb and 23 Na have interaction 
parameters within this range. At 8 = 1, {g\ = gi = g\2 — g), the aspect 
ratio is \/3- For indistiguishable components, the combined lattice would 
then appear triangular. 

As described in Sec. 3, the condition of phase separation for nonrotating con- 
densates is given by Eq. (29) (i.e., S > 1). The presence of a vortex lattice naturally 
modulates the density of each component, with the high density regions of one 
component overlapping with the low density regions of the other. Thus, the system 
is effectively phase-separated whenever staggered vortex lattices are present, even 
for 6 < 1. In particular, the vortex lattice near 6—1 (Fig. 16(e)) is made up of 
alternating rows of vortices of each component, and the system therefore contains 
stripes in which one component has a high density and the other component has a 
very low density. 
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5.2. Numerical simulations for rotating two-component BECs in 
the Thomas-Fermi regime 

This section addresses a full numerical study of the GP equations on the vortex 
lattice structure in two-component BECs 88 . The assumptions in the LLL calculation 
shown in the last subsection are that (i) the rotation frequency is assumed to be very 
close to the radial trapping frequency and that (ii) a vortex lattice is assumed to be 
perfectly regular 110 . Therefore, a vast experimentally accessible parameter regime 
of the lattice configuration remain to be explored. The full numerical simulations of 
Eq. (46) and (47) reveal a rich variety of vortex structures which have eluded the 
analytic treatment. Here, we show some detailed results of our simulations. 

We numerically calculated the equilibrium solutions of Eq. (46) and (47) through 
the conjugate gradient method 23 ' 172 . Following Mueller and Ho, we also assume an 
identical trapping potential for both components (V\ = V2 = V) and the identical 
intracomponent interactions g\ = g^. We confine ourselves to the 2-D problem. The 
dimensionless coupling constants (defined in Eqs. (46) and (47)) are u\ = 1*2 = 
u = 4000 under the assumption N± — N2; the normalization of each wave function 
is taken as J d?r\ipi\ 2 = 1/2. The intercomponent interaction U12 and the rotation 
frequency fl are treated as variable parameters. Figure 17 shows the numerically 
obtained phase diagram of vortex states in the (8 = Ui2/u,fl) plane. The upper 
limit of the rotation frequency is SI — 1 , because for il > 1 the centrifugal potential 
dominates the harmonic trapping potential and therefore a BEC cannot be trapped. 
Also, we do not discuss the slow rotation regime below 51 = 0.40. The criterion of 
phase separation 8 = 1 is reflected on the structure of vortex states, as explained 
below. 

In the overlapping region d < 1, two types of regular vortex lattices are obtained 
as the equilibrium state. For 8 = 0, the two components do not interact; therefore, 
triangular vortex lattices form as that in a single component BEC. As 8 increases, 
the positions of vortex cores in one component gradually shift from those of the 
other component and the triangular lattices become distorted. Eventually, the vor- 
tices in each component form a square lattice rather than a triangular one, which is 
consistent with the result of the LLL calculation done by Mueller and Ho 87 . How- 
ever, Fig. 17 shows that the stable region of the square lattice depends not only on 
8 but also on fi. Figure 18 shows the structural transition in closer detail. While an 
increase in 8 indeed causes the deformation of the lattices from triangular to square, 
the transition occurs at a significantly higher value of 8 than that of the LLL result 
(8 = 0.373). For example, for il = 0.7 the transition occurs at 8 ~ 0.65. This implies 
that an increase of rotation frequency also causes the transition from triangular to 
square lattices. From Fig. 18, the lattice changes from triangular to square around 
O = 0.7 — 0.8 for 8 = 0.6. We find that the two vortex lattices arc interlaced in such 
a manner that a peak in the density of one component is located at the density 
minimum of the other. As a result, the total density = \4>i\ 2 + \ip2\ 2 obeys the 
Thomas-Fermi distribution applied to the overlapping two-component BECs with 
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Fig. 17. Q-S phase diagram for vortex states in rotating two-component BECs. (a) Table of data 
for the phase diagram in (b). A: triangular lattice, O: square lattice, eg): stripe or double-core 
vortex lattice, o: vortex sheet. Because of the continuous change between triangular and square 
lattices, their boundary is shown by a combination of A and o. The plots at Q = 1 show the results 
obtained by Mueller and Ho [Phys. Rev. Lett. 88, 180403 (2002)] based on the LLL approximation, 
(b) Vortex phase diagram. The boundaries of the vortex sheet phase arc given by the condition 
b = A p (solid curve) and 2b = -Rtf (dashed curve) (see text). 



rigid body rotation pt(t) = y/2a/n — ar 2 with a = (1 — fl 2 )/u(l + 5). We have 
confirmed that this generic feature holds for other parameter regimes. 

The pseudospin and the NLctM provide insights for understanding why the 
square lattice is stabilized in two-component BECs. According to the energy func- 
tional of Eq. (14), two components interact via the intercomponent interaction 
<7i2|V'i| 2 |V'2| 2 ; however, the velocity field in one component is independent of that of 
the other. Therefore, the intercomponent relation is determined only by the density 
distribution of the condensates that minimizes the interaction energy. According to 
the NLcrM (48) , we can rewrite the interaction energy in terms of the total density 
Pt and the z-component of pseudospin S z = |xi| 2 — | X2 1 2 - In this representation, the 
spin-up component of S z corresponds to the density peak of ipi at the vortex core of 
ip2, and vice versa for spin-down components. Now, we have assumed wi = u 2 = u, 
so that ci = and the interaction energy reads 



2 Up T 

ar — - 

4 



' ■(1 + S) + (1-6)S1 



(114) 



This expression is similar to that of a binary alloy, in which N/2 "A" atoms and 
N/2 "B" atoms are distributed on TV lattice sites 88 . A larger 5 makes a smoother 
total density px that is more favorable for reducing the first term of Eq. (114). 
This results in a shift of the positions of vortex cores of each component. Therefore, 
the spin-spin interaction c 2 in the second term plays a crucial role in stabilizing 
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Fig. 18. Two-dimensional density profile of the condensates ipi and tp2 with m = 4000 in the 
overlapping region 5 < 1. We choose typical values of rotation frequency Q = 0.6, 0.7, 0.8, and 
those of the intercomponcnt interaction strength S = 0.4, 0.6, 0.8. The plot range of the density 
profile is [—W<x,y< +10] in units of fer- 



tile square lattice. When the coefficient 1 — 5 is positive, ci > and the interac- 
tion between spins is anti-ferromagnetic, which favors a square lattice because the 
triangular lattice can become frustrated. This can be understood by representing 
the vortex lattice state with the pseudospin. The vortex lattices of two-component 
BECs can be regarded as a square lattice of skyrmions, in which skyrmions with a 
radial-disgyration texture (+z) alternate with those having a hyperbolic-disgyration 
texture (— z). The profile of S z in Fig. 19(b) shows a square lattice structure, which 
is equivalent to the "continuous" Ising model of spins; actually, the second term of 
Eq. (114) is equal to the first term of the Ising Hamiltonian in a continuous limit 

#s P in = J Hi.j S z {vi)[S ' z {ri) + (r j - Vi) • V 'S 'z(rj) H ]. As both S and tt increase, the 

interlocking of vortex lattices becomes stronger and the anti-ferromagnetic nature 
is made more pronounced, which lead to a nontrivial boundary between triangular 
lattices and square lattices in the phase diagram of Fig. 17. 

As S exceeds unity, ci becomes negative and the system enters a ferromagnetic 
phase. Then, the condensates undergo phase separation to spontaneously form do- 
mains having the same spin component. Concurrently, vortices begin to overlap at 
5 <~ 1. For 6 = 1, an SU(2) symmetric case, interesting structures appear as shown 
in Fig. 20. In Fig. 20(a), each condensate density forms a stripe pattern such that 
the lines of vortices for component one overlap those for component two. Mueller 
and Ho analytically derived the lattice structure in Fig. 20(a), with the assumption 
of a perfect lattice 110 . However, our numerical treatment shows that the lattice is 
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Fig. 19. Spin texture for the square lattice solution 8 = 0.8 and Q = 0.7 (cf. Fig. 18). (a) The 
spin texture profile projected onto the x-y plane. The positions of the defects with +z and — z arc 
marked by Q and ®, respectively, (b) The spatial distribution of S z . 




y y 



Fig. 20. The density profiles of the condensates ipi (left) and i/>2 (right) for Q = 0.7 and 6 = 1, 
i.e., u\ = «2 = "12- The profiles (a) and (b) were obtained by numerical simulation starting from 
different trial wave functions. The energy deference between these states is AE ~ 10 — 5 fiiij±N. 

significantly different when the energy is only slightly different [Fig. 20(b)]. For the 
same parameters, we can obtain another equilibrium state called a "double-core 
vortex lattice" in Fig. 17, where a vortex lattice of component 2 is made up of pairs 
of vortices with the same circulation and the vortices in component 1 surround those 
pairs. Therefore, various metastable structures will be observed in this parameter 
region. 

For strongly phase-separated region 6 > 1 (c2 < 0), the domains of the same 
spin component, at which the other-component vortices are located, merge further, 
resulting in the formation of "serpentine" vortex sheets. A typical example is shown 
in Fig. 21(a). Singly quantized vortices line up in sheets, and the sheets of component 
1 and 2 are interwoven. Figure 21(b) shows that each superfluid velocity (i — 
1,2) jumps at the vortex sheet, following approximately the velocity of solid-body 
rotation. In the region 5 > 1, even though the value of S is changed, the total density 
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Fig. 21. Vortex sheets, (a) The density profiles of the condensates for S7 = 0.75 and <5 = 1.1. The 
vortex sheets arc made up of chains of singly quantized vortices whose positions are marked by X . 
(b) Bold curves show the density profiles of \J>i and ^2 hi the radial component, and thin curves 
the corresponding velocity profiles |vj| = |Vargi/>,|. The velocity of rigid body rotation v s = Qr is 
shown by a dotted line for comparison, (c) Model of the vortex sheet state. See the text for details. 

is fixed by the Thomas-Fermi distribution with 5 = 1. 

A stationary vortex sheet has been observed in rotating supcrfhiid 3 Hc-A 173 . In 
that system, the vortices are bound to a topologically stable, domain-wall soliton 
across which the direction of the orbital angular momentum of the Cooper-pairs 
have opposite orientations. Moreover, experiments with rotating supcrfluid 3 He- 
A show that the equilibrium separation between the sheets is determined by the 
competition between the surface tension a of the domain wall and the kinetic energy 
of the superflow 173 . 

How do these properties differ for two-component BECs? To estimate the sheet 
separation in two-component BECs, we consider a simple model shown in Fig. 21(c). 
In this model, each velocity is assumed to be uniform between the vortex sheets of 
the same component; specifically, the value of Vi increases by 206 across every sheet. 
The total density px = \ipi\ + \ip2\ is constant, and the domain boundary with 
the penetration depth A is approximated by the linear profile. We calculate the free 
energy F = E — HL Z in the range < r < 2b with the sheet distance b being deter- 
mined so as to minimize F per unit area. Although the condensate density is lower 
near the edge because of the trapping potential, this model should be valid in the 
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central region. We first estimate the optimized penetration depth A p by minimizing 
the surface tension a of a single domain wall with respect to A 136 , which is the sum 
of the quantum pressure energy E qp /2Trb w (px/2A) and the interaction energy 
in the overlapping region £'i nt /27rfe = u(S — l)h.p\/12. Minimization of the surface 
tension with respect to A yields A p = y/6/u(5 — ljpT- Then, the corresponding 
surface tension is written as 



3/2 U(5- 1) 

o" = Pt Y g • ( 115 ) 



Next, the flow energy per area in rotating frame is 

r 2 & /■„,. _ 0^2 O0„_02?,2 



1 2 ^ - »r) 2 29 /OT ^& 2 

i^x dr 4> — 2 — =^6^' (n6) 

where b 3> A p is assumed. Thus, the free energy per unit area is 

F 29p T Q 2 b 2 2a up 2 
4tt6 2 768 b 2 ' { 1 

Minimizing this energy with respect to b, one obtains 

, / 768,7 \ 1/3 / 768 \ 1/3 

By using the Thomas-Fermi density at r — with 5 = 1 as the value of pt in 
the denominator of Eq. (118), the sheet spacing b oc (5 - - Q 2 ) 1 / 12 /^ 2 / 3 

is consistent with what we found numerically; for example, for parameters used in 
Fig. 21 we obtain b — 3.096ho- The vortex sheet is expected in the region b > A p 
and 2b < i?TF = as shown in the shaded region in Fig. 17(b). When 

2b > i?TF (i-c, when b becomes comparable with the size of the condensate), a 
strong ferromagnetic interaction reduces the size of the domain wall region and the 
serpentine structure disappears as shown in Fig. 22, realizing rotating "droplets" in 
which most vortices are located in the low density region. 

What structures can be observed experimentally? The two-component BECs 
realized in JILA are a mixture of the states |1, —1) and \2, 1) of 87 Rb. This mixture 
has scattering lengths in the ratio of an : 022 : a.12 = 1 : 0.94 : 0.97 and so 5 <~ 1 for 
equal particle numbers. Therefore, lattices with partially overlapping vortices, as 
shown in Fig. 20, are expected to be observed. 85 Rb atoms in the |2, 2) state make 
an interesting system because both the intcrcomponent scattering length and the 
scattering length of 85 Rb |2,2) can be controlled by using the Feshbach resonance. 
The MIT group has observed the phase separation of a mixture of 23 Na BECs with 
1 1, 0) and |1, 1) state 72 , which has S > 1 (ai : a 2 : a 12 = 1 : 1.035 : 1.035). A mixture 
of 41 K and 87 Rb BECs which is reported recently 23 lies deeply in a phase-separated 
region. The vortex sheets should therefore be observed at high rotation frequencies. 
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Fig. 22. Density profiles of the condensates for f2 = 0.55 and <5 = 1.5. Positions of vortices are 
marked by X . 

5.3. Effect of internal coherent coupling on lattice structure 

As discussed in Sec. 4.3, an external coherent-coupling field makes bound vortex 
pairs in a two-component system, inducing an effective attractive interaction be- 
tween them. A question then naturally arises: what happens when an external field 
is applied to rapidly a rotating two-component BEC? Since two components inter- 
act with each other through their relative phase, the structure of the vortex states 
should be greatly affected by this effect. Zhai et al. studied a similar situation that 
involves the structure of a vortex lattice in a condensate trapped in a double-layered 
potential 174 . The interlayer coherence, which seems to be a small contribution in 
such a situation, was instead shown to have a significant effect on the structural 
transition of the lattices. Their treatment was in the limit U12 <C u, and thus their 
analyses were simpler than the case discussed in this subsection. Here we present 
some numerical results about the structure of vortex lattices with a coupling field. 

First, we discuss the simple case u\ = u 2 = u i2 {c\ = c 2 = 0) with the SU(2) 
symmetry. As the rotation frequency increases, a large number of vortices nucleate 
and form stripes or double-core lattices as shown in Fig. 20. When wr is turned on, 
the vortices begin to form pairs and then form a lattice of meron pairs. The resulting 
equilibrium structures for f2 = 0.8 and several values of wr (0, 0.2, 0.5) are shown 
in Fig. 23(a). The meron pairs are characterized by their axisymmctric topological 
charge and vorticity, yielding a their square lattice. This lattice of meron pairs 
is equivalently a lattice of skyrmions under the basis transformation (ipi,ip2) — > 
(ip+,il>-) as described in Eq. (58) and shown in Fig. 23(b). As wr increases further, 
vortices of one component overlap completely with those of the other, forming a 
triangular lattice as in a single-component BEC. Viewed in terms of the ip± -basis, 
the decrease of the particle number in the ip- component effectively weakens the 
mutual repulsion between the two components, thus giving rise to a transition from 
square lattices to triangular lattices. Therefore, application of an external coupling 
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Fig. 23. (a) The density profiles of the equilibrium solutions ipi and ip2 for Q = 0.8 and loh = 0, 
0.2, 0.5 from top to bottom, (b) The corresponding profiles of i/>+ and obtained by the basis 
transformation t/>± = (i/ji ± tp2)/V2. 

drive could be an excellent tool to explore a rich variety of vortex structures in 
two-component BECs. 

However, in the case that u\ ^ U2 ^ u\2, the meron pair is characterized by an 
anisotropic topological charge as shown in Fig. 12. A consequence of this anisotropy 
is that the resultant vortex state exhibits a distorted lattice structure. Typical ex- 
amples are shown in Fig. 24, where (a) and (b) correspond to the antiferromagnetic 
and ferromagnetic cases, respectively. From the bottom panels in Fig. 24, the direc- 
tion of the ordering depends strongly on the distribution of the topological charge in 
each meron pair. The collective mode of this distorted lattice is expected to exhibit 
anisotropic wave propagation, depending on the polarization of the meron pairs. 

5.4. Experimental observation of square lattices in two-component 
BECs 

Schweikhard et al. observed well-separated vortices forming in two-component BECs 
that later formed ordered square lattices 165 . As in Matthew et aFs experiment 11 
mentioned in Sec. 4.1, the two-component system consists of two hyperfine levels 
of the 87 Rb atom: \F = l,m F = -1) = |1) and \F = 2,m F = 1) = |2>. To study 
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Fig. 24. Vortex lattices of SU(2)-symmetric BECs with internal coherent coupling. The upper 
panels are the profiles of the condensate density \ipi\ 2 and |V>2| 2 and the lower, larger images are 
the relative phases for SI = 0.80 and uj-p. = 0.2, the centers of mass of the molecules are linked by 
dashed lines. The bottom panels show the distributions of the topological charge <?(r). (a) C2 = 20 
(antiferromagnetic) . (b) C2 = — 20 (ferromagnetic). 

the rotational properties, they created a regular triangular vortex lattice in a BEC 
in the state |1) that contained (3.5 — 4) x 10 6 atoms and rotated at a frequency 
VL w 0.75 x u)± about the z-axis. Then, a short pulse of the coupling drive transfered 
a 80 — 85% fraction of the population into state |2). After a variable wait time, they 
took two images of the condensates. One was a nondestructive phase-contrast image 
[Fig. 25(d)], which was taken along an axis perpendicular to the rotation axis. The 
other was a destructive image to resolve the vortex core structures [Fig. 25(a) and 
(b)] taken along the rotation axis after the free expansion. 

Figures 25 (a) and (b) show the formation and decay dynamics of the square lat- 
tice structure of the expanded |1) and |2) states. For the first period 0.1— 0.25 s, there 
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Fig. 25. Experimental observation of square lattices in two-component BECs. (a) Time sequence 
of images of state |1), after ~ 80% population transfer to |2), showing the evolution from a 
hexagonal lattice over a turbulent stage to a square structure and back to a hexagonal lattice, 
(b) Similar to (a) except for state |2). In this case, a square lattice forms (iv) before the lattice 
decays, (c) Time evolution in reciprocal space. Intensity in an annulus along the <p coordinate. 
The ordinate is the angle. The initial 6-peak structure of the hexagonal lattice vanishes quickly 
due to turbulence. The four peaks in 3 — 5.5 s indicate a square lattice. During 6 — 9 s a transition 
back to hexagonal lattices occurs, (d) Two-color edge-on images of the initial turbulent evolution. 
State |1) (|2)) appears bright (dark) on gray background. The fine filament structures in (ii-v) are 
due to mutual filling of vortex cores. [Schweikhard et al, Phys. Rev. Lett. 93 210403-2 (2004), 
reproduced with permission. Copyright (2004) by the American Physical Society.] 

is little dynamical behavoir in either component and certainly no structural transi- 
tion in the vortex lattice. From 0.25 — 2 s, turbulent behavior appears in both com- 
ponents in which vortex visibility degrades significantly, as shown in Figs. 25[a(ii)] 
and 25[b(ii)]. This turbulence is a direct consequence of the transition from over- 
lapping hexagonal vortex lattices to interlaced square lattices. From 2 — 3 s, square 
lattices emerge from the turbulent state (Fig. 25[a(iii)]). From 3 — 5.5s stable square 
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lattices are observed in both components (Figs. 25[a(iv)] and [b(iv)]). At this stage, 
around 4 s, despite the large (80 — 85 %) initial population transfer to the state |2), 
the number of |2) atoms has decreased to only 1.5 x 10 5 , while the state |1) contains 
(5 — 7) x 10 5 atoms. This is due to a larger inelastic scattering loss of the |2) state 
than that of the |1) state. As the |2) state population continues to decay, the vortex 
lattice planes bend (Fig. 25[a(v)]) and a transition back to a hexagonal lattice in 
the state |1) takes place (Fig. 25[a(vi)]). No turbulence occurs during the transition 
from a square to a hexagonal lattice. 

To discuss the structural transition more clearly, Fig. 25(c) shows the time evo- 
lution of the intensity of the reciprocal lattice peaks of the vortex lattice. In this 
plot, the triangular lattice has six peaks at <j) = 60° x n with integer n, whereas the 
square lattice has four peaks at 90° x n. Initially six peaks are visible, which corre- 
sponds to the reciprocal lattice of a triangular vortex lattice. Because of subsequent 
turbulence these peaks vanish between 200 ms and 2 s. After 2 — 3 s a four-peak 
structure appears, which is a clear signature of a square lattice. At 5.5 s two addi- 
tional peaks appear, which signals the onset of the transition back to the triangular 
lattice. This experiment thus demonstrated that, although the state is transient, 
the intercomponent interaction yields the square lattices of vortices as predicted in 
Refs. 87 and 88. 

Figure 25(d) shows a time sequence of side-view images, where the |1) (|2)) state 
appears bright (dark). As can be seen in Fig. 25[d(i)], a spatial phase separation 
from the initially overlapping two components to a ball-shell structure occurs within 
50 — 100ms. During the period of this axial separation, both the individual vortices 
and the overall vortex lattice remain remarkably stable, as viewed along the rotation 
axis (see Fig. 25[a(i)] and Fig. 25[b(i)]). Around 600ms, fine filament structures 
appear at the intercomponent boundary (Fig. 25[d(ii)]), concurrently with the full 
development of vortex turbulence seen in Fig. 25[a(ii)]. The filament structures 
distort and fill out the whole BEC as the vortex turbulence peaks at around 1 s 
(Fig. 25[d(iii)]), and straightens up at around 2 s (Fig. 25[d(iv)]), concurrently with 
the restoration of vortex visibility in top view images. Subsequently, the filaments 
become less visible as the state |2) decays (Fig. 25[d(v)]). These structures are 
interpreted as vortex cores in either component being filled by fluid of the other 
component, forming a skyrmion lattice. In the experiment, interference techniques 
were used to investigate the detailed dynamics of the phase separation. Also, the 
stability of square-lattice structure was confirmed by observing the fast relaxation 
of the Tkachenko excitations, as in the single component BEC 42 ' 45-49 . 

6. Conclusions 

We have reviewed the physics of quantized vortices in rotating multicomponcnt 
BECs, with an emphasis on unconventional vortex states in two-component system. 
The vector nature of the order parameter gives rise to intriguing structures and 
dynamics that are absent in a single component system. Because analogous physical 
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phenomena and theoretical formulation are found in other physical systems which 
involve multicomponent order parameters, an indepth understanding of quantized 
vortices in multicomponent BECs may offer physical insight in other systems. 

For the slowly rotating regimes, we discussed the structure of single vortex states 
in two-component BECs with and without internal coherent coupling by using both 
numerical simulations of the coupled CP equations and the variational analysis 
based on the nonlinear a model. The vortex structure is significantly altered from 
that of a single-component BEC because of the intercomponent interaction. We 
showed that a pseudospin is a useful concept to represent the vortex state, which 
reveals structures and energetics of the important topological excitations called 
skyrmions and merons. A vortex molecule (a meron pair) is predicted to be stabilized 
in the coherently coupled two-component BEC in a rotating potential, and the con- 
tributions that breaks SU(2) symmetry make the vortex molecule have anisotropic 
vorticity. We showed the theoretical studies on the 3-D nontrivial skyrmion, which 
is regarded as a composite vortex ring and is closely related to a proposed topo- 
logical excitation in the early universe. Also, we gave an introdutory discussion 
on quantized vortices and monopolcs in spin-1 BECs. Much remains to be inves- 
tigated concerning the mutual interactions between exotic topological defects and 
their nontrivial dynamics. 

We described how a fast rotation can stabilize a lattice of skyrmions and merons. 
Such a lattice consists of vortices in each component. The analytical and numerical 
studies revealed a rich variety of vortex phase diagrams in the variables of rotation 
frequency (0.4 < fl/w± < 1) and an intercomponent interaction strength. The inter- 
component repulsive interaction causes one lattice to displace relative to the other, 
stabilizing mainly the interlaced square lattices in an antifcrromagnctic regime and 
the serpentine vortex sheets in a ferromagnetic regime. We proposed that an exter- 
nal driving field is a useful tool to control the strength of the interaction between 
the two components and, thereby, provides a unique means to study the rich variety 
of vortex phases. Much work is necessary to reveal further details on the lattices 
of skyrmions and merons in rapidly rotating two-component BECs; moreover there 
remain further interesting problems such as structural transitions and collective os- 
cillations of the lattice of skyrmions and merons, and new phenomena related to 
double-layer quantum Hall physics 4 . 

Another type of two-component system of ultra-cold atoms is boson-fermion 
mixtures. Such studies, which includes 7 Li- 6 Li (Refs. 175, 176), 23 Na- 6 Li (Ref. 177), 
and 87 Rb- 40 K (Refs. 178, 179, 180), are analogous to the system of 3 He- 4 He interpen- 
etrating supcrfluids 181 . In a spin-polarized fermionic system, s-wave fermion-fermion 
collisions are forbidded by the Pauli exclusion principle, and the Fermi pressure is 
the dominant contribution to the ground state properties. The boson-fermion mix- 
ture will provide another field to investigate rich vortex phenomena mediated by 
the boson-fermion interaction, which is repulsive 175 ' 176 or attractive 178 ' 179,180 , and 
can be tuned through the Fcshbach resonance 182,183 . The stability of vortex states 
in a 87 Rb- 40 K mixture was studied theoretically in Ref. 184. Quatized vortices in 
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spin-2 BECs are also an intriguing subject and will soon open up a new research 
field because spin-2 BECs have been successfully loaded in an optical trap 74 ~ 76 . 
This system is even more attractive than the spin-1 system, because it is expected 
to have a rich variety of novel ground state phases, a magnetic response, and spin 
dynamics 185-187 . A method of topological phase imprinting 13 ' 14 ' 112 is also a power- 
ful tool to generate a quantized vortex 188 and the experiment toward this direction 
is in progress 189 . 
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